{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 543,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 544,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAENCAYAAAD0eSVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF+JJREFUeJzt3X2UZHV95/H3VxBK1IgsojyIQslKjBLxoC5HwUpwp6Kd\nEHCDGqJskLgxq9XtEDUEIUzMCp51Ccl0ohCCrg8bPPiEIYVtGx0wiRsXzYCIaLB1VIiMD0GRQI2B\n+e4f99ZMTVvTT9M9t/r2+3VOnbr3V1X3fu+dmc/86verqhuZiSSpXh5WdQGSpOVnuEtSDRnuklRD\nhrsk1ZDhLkk1ZLhLUg0Z7pJUQ4a7JNWQ4S5pWUXEUUPaDo2IA6qoZ60y3KU9EBHby9uRVdcyCiLi\naOA/DXnoe8Cb9nI5a5rhXmNR2DIQQMdWXdNKiYijI+IjEfHdiHggIr4dER8vw2Yl/SnwJ8CPV2Lj\ns/78HoyI70TE+yPikCVs4+SVqHGW387Mq2c3ZuaDQDciztoLNQjDve5OAgZ7lK9cyZ1FxH4ruf15\nfBQ4DbgFeBewmaIH+YSlbGyhx5KZ6zPz3My8Zyn7WYTrgMuBBM6k+A9lMRKI5S5qUET8PHDnrLYz\nI+J7AJl5E/DClaxBOxnu9faK8n5zeX9m/4GIeE/Zmzt3oO3dg20RcWREfCAi7oqIeyLiExHxcwPP\n7/coXx8R3wBuL9v/qnxNLyLujYhPRcTTB173/Ii4NSLui4j3lvvYHhGXLWS/s0XEQcAzgHsy8z9n\n5msz81TgkP6xL+VYIuJ/L+Ac7TIsExFHlOf2m+XxfzkiTljKcc1yVWa+DjivXD9uoKbdnu+I2ELx\nH3wAm8pa/+sCz0n/+C9bYI2/DHx6VtvtwA0D69+LiKcscHvaE5nprYY3YH/gX4GHgFa5vB04uXz8\nF8v1fyzXHw7cA/wEeDxwAPC18vVd4D3ANmArcFD5mu3l7d+AdwPvKNv/Hngf8GcU/7C3A18uHzuw\n3M924FPAFPBguf7HwCPm2O9/2M2xPhy4t9zGzeV2fhV4ZPn4ko5lvnM08LqHKAL0AOCfy7bbKXra\nNwK/Mk8NQ4+r3P6Wcnu/CuwHbCzX/2rgOXOd7wuBH5Vt15Tn5oSF1FOeh4eAP17g37lrgZjV9rvA\n7wysnwW8rOp/H2vhVnkB3lboDxZeUv6D/g5Fr+095fpflI8H8K3yH++TgbHy8evLx88o178NXEYx\nDPC1su23y+f0A/E3Z+37MKADXAJMDgTgoRTvJrYDdww8fzM7w33e/e7meM9g538a/dt3yiBb0rHM\nd44GXtcP918r1+8EGgPP2WcPjmvLrGPaDmyi/E9pnvP9hFnbOHnW+ZrvnDwB+I/A42Ydy98PrF8F\nNMvl6SH1Xwc8dWD9V4Dfrfrfx1q47Yvqqj8kc11mZkRcSzHmfkZEvC4zfxIR7wN+H3gp0H87/p7y\n/snl/eHARLmc5a05sJ8E/qG/EhHHAP8EPHJITY8rtwfw1YH2rwA/TxGmT1rgfneRmR+MiI8BLwBO\nBl5N8Q7kQuDvlnIs5Xmb6xwNCqD/EcBbM7M3sJ2HIuLJSzmuAddRvMM5HXgucCzw2TnOd1Cc77t3\ns71568nMu4e8/kTgm1BM2AMnZuZM+dg+uxQQsQ/wlMwc/LN+gOIdiFaYY+41FBEHAi8uV38rIrYD\nHy7XH0PRewJ4b3n/G8CpFG/fry3bvlHefz4zH1be9gEOAi6etcttA8tjFEGzudzX7AnN/oTbMQNt\n/U/xJEUvc6H77R/vvhHx/Mz8SWZ+MjMvBN5WPvyoRW5z26z1uc7RoAS+Xi4fFxGNgfr2YXHnc5ir\nMvO/UAyVNCiGZ2D353twAvWh8n4wfOetJ4rPph8bEY8beN0vAZ8ol48Hbh147MFZNT8buGlW22Mo\nhgi1wgz3enopRe/oXoog6t/uKB9/JUDZo7qJYjLyMcCHMrMfbtdTBMAJEfHZiLg8Iq4H/oWByTx+\n+hMY/Z7eUymGCKZnPe9vgB8Cx0TE30bEVLn/voXud1AD+ExE3BYR/ycirgQuKB/75B4cy3znaPbr\nrqc4x4cCmyPiioi4kSKAl3Jcw/whRVg/KyLWMf/5hmJoCeAtEXFZRByxwHouAb4MnD+wrTbFOwXK\n4/pURJxart8dEY8aeO5JwA0RcfpA26EUwz9aYYZ7Pf0GRc/t8sx8Sf9GMVQB8EsR8dhyuT/EkOzs\npZKZ91NMKF4NPJFiIuwYiom7wbfZs6/TeA3FOOxD5evfOvi8zPwRxTuHL1EML2wF/rp8zrZF7HfQ\nAxTjxtuAF1EMSf0r8Bbgf+7BsfQNPUezX5eZDwCnlNs9oNzPwcBdSzyu/rZ31JWZ3yprSOD3mOd8\nlzZQBOqJwDhwyALryYEbZQ/+SODUiHgxcD/F0M/95fNvBJ4zsN9/Bo5g1576MxkY+tLKiUyvoaq9\nKyJ+JjPvLZcfBtxG0fM8JzPfXWlx2q2IeAXws5n55t08fiDwhsy8YDePN4CLM/PcYY9reTmhqipc\nFREPUnxc8AUUwX4XO+cFNJqew/DJZAAy84cR8f2IODgzvz/kKS8Hrlix6rQLe+7a6yLiPOC/U3zJ\naCvF2/QLMvPrc75QI6/8BM1vZeaVs9qfCDwrMz9WTWVrj+EuSTXkhKok1dDIjLlHRPouQpIWbegP\nwtlzl6QaMtwlqYYMd0mqIcNdkmrIcJekGjLcJamGDHdJqiHDXZJqyHCXpBoy3CWphgx3Saohw12S\nashwl6QaMtwlqYYMd0mqIcNd0qrR7XZpt9u0Wi3a7TbdbrfqkkbWyFysQ5Lm0u12mZiYYGZmZkdb\nf3lsbKyqskbWyFxD1SsxSZpLu91menp6aPvU1FQFFY2MoVdisucuacEihuZIpXq9XtUljCTH3CUt\nWGZWdlu3bt3QmhqNxl4+C6uD4S5pVRgfH6fZbO7S1mw26XQ6FVU02hyWkbQq9CdNJycn6fV6NBoN\nOp2Ok6m7sewTqhGxf2ZuG9LeyMzdDo45oSpJSzJ0ImRZh2Ui4peBR0fEqRHxmVkPHxERL1zO/UmS\nhlu2cI+IQ4GfyczvA3cA/3fw8cz8GvC0iHjEcu1TkjTccvbczwY+Wi6fCHx2yHO6wK8v4z4lSUMs\nOtwj4sUR8cpy+a0RcWT50CGZ+UC5/Bxgc0S8JCK+0H9tZs4Az9jToiVJc1tKz/0UYHO5/KzM/Fa5\nPPhh06cBz87MjwDPn/V6P6EjSStsKUH7jMz8UkTsD/xkoP3hABHxqHL99IjYnpkfnfX6A3a34Q0b\nNuxYbrVatFqtJZQnSVpUuEfEAewM5+cCN0fEyZn5GeChsv3ZwN8AU8AZEbEtM68f2Mz23W1/MNwl\nSUu32J77c4HHRMQYcBCwP/Dv5WP3l/fHApuAO4FHAD/qvziKH6b48Z4ULEma32LD/XnA6zLzxiGP\n3RkRj83Mdw60vWHWc44DPrfIfUqSFmmxE6pHM+vz6wOuBM6Y5/WnAB9c5D4lSYu0rD8/EBEnAd8c\n+ATN4GM/B+ybmbfs5rX+/IAkLd7Qnx/wYh2StLqt/G/LSJJGg+EuSTVkuEtSDRnuklRDhrsk1ZDh\nLkk1ZLhLUg0Z7pJUQ4a7JNWQ4S5pXt1ul3a7TavVot1u0+12qy5J8/CqSJLm1O12mZiYYGZmZkdb\nf3lsbKyqsjQPf1tG0pza7TbT09ND26empiqoSLMM/W0Ze+7SKlFc62Z09Hq9qkvQHBxzl1aJzKzk\ntm7duqH1NBqNvXwGtBiGu6Q5jY+P02w2d2lrNpt0Op2KKtJCOCwjaU79SdPJyUl6vR6NRoNOp+Nk\n6ohzQlWSVjcv1iFJa4XhLkk1ZLhLUg0Z7pJUQ4a7JNWQ4S5JNWS4S1INGe6SVEOGuyTVkOEuSTVk\nuEtSDRnuklRDhrsk1ZDhLkk1ZLhLUg0Z7pJUQ4a7JNWQ4S5JNWS4S1INGe6SVENzhntEbIyIrRGx\nPSKuG2jfUrb1b5vL9sMj4qaIuDciLh14/vqIuGLlDkOSNGi+nnsCVw8sD7bfCLy8vL2pbH8tcBBw\nKbA+IpoRcTAwDpy/XEVLa0m326XdbtNqtWi323S73apL0iqw71wPZuZERDyJIpwHBbAFuD4z7xto\nPwC4G9gEXAQ8GngDsDEzf7BcRUtrRbfbZWJigpmZmR1t/eWxsbGqytIqsJAx9xjSlsBZwL3lsM2r\nyvargeOBG4CbyteeDGzc81KltWfjxo27BDsU4T45OVlRRVot5uy5z+FK4CvA/sDbgCsi4tOZ+bmI\nOAo4ErgF+DjwRuA1EXEu8D3g7My8fdhGN2zYsGO51WrRarWWWJ60fCKG9W+q1ev1qi5BI25J4Z6Z\nF/eXI+IEYD1wDLAlM7cCWyPiNGAbRS/+Y0ALOAe4EDhz2HYHw10aFZk5/5NWSLvdZnp6+qfaG41G\nBdVoNZkz3CNiDHh6uXpkRJwDfB74HxS98n0phmfuB24deN1+wMXA6cA+FMMzLwOeCdyxvIcg1df4\n+DgzMzO7DM00m006nU6FVWk1iLl6JRGxCXgBxRh7lPdvBk4CnkMxgXob8ObM/OTA694IHJaZ68v1\n84DzKCZbX56ZNw/ZV1bZQ5JGVbfbZXJykl6vR6PRoNPpOJmqQUPHDecM973JcJekJRka7n5DVZJq\nyHCXpBoy3CWphgx3Saohw12Sashwl6QaMtwlqYYMd0mqIcNdkmrIcJekGjLcJamGDHdJqiHDXZJq\nyHCXpBoy3CWphgx3Saohw12Sashwl6QaMtwlqYYMd0mqIcNdkmrIcJd2o9vt0m63abVatNttut1u\n1SVJC7Zv1QVIo6jb7TIxMcHMzMyOtv7y2NhYVWVJCxaZWXUNAEREjkotUrvdZnp6emj71NRUBRVJ\nuxXDGu25a6RFDP17W5ler1d1CdKCOOaukZaZldzWrVs3tJ5Go7GXz4C0NIa7NMT4+DjNZnOXtmaz\nSafTqagiaXEclpGG6E+aTk5O0uv1aDQadDodJ1O1ajihKkmr29CJKYdlJKmGDHdJqiHDXZJqyHCX\npBoy3CWphgx3Saohw12Sashwl6QaMtwlqYYMd0mqoTnDPSI2RsTWiNgeEdcNtD8vIr4YEb2I+EJE\nHF+2Hx4RN0XEvRFx6cDz10fEFSt3GJKkQfP13BO4emCZiGgAHwYeCbweeDzwoYh4GPBa4CDgUmB9\nRDQj4mBgHDh/+cuXJA0zZ7hn5gRw2azmFwGHAO/IzMuBdwFHAS3gAOBuYFP53EcDfwRszMwfLF/Z\nkqS5LGTMffYvjh1V3t9V3t850H41cDxwA3BT+dqTgY17VKUkaVGWY0J1R/hn5ucoQv65wEnA/wLe\nCLwmImYi4h8j4meXYZ+SpDks5WId3yjvn1jeH17efx0gM7cCWyPiNGAbRS/+YxTDNucAFwJnDtvw\nhg0bdiy3Wi1ardYSypMkzXmxjogYA54OXAJ8EZgE/h/wSeB+4O3ABUAPeEr/ahsRsR9wM3A68C/A\nPcA7gOcDd2Tmy4bsy4t1SNLiDb1Yx3zhvgl4AcUnZaK8P5ui9/7nwFOBLwGvzsx/GnjdG4HDMnN9\nuX4ecB7FZOvLM/PmIfsy3CVp8RYf7nuT4S5JS+Jl9iRprTDcJamGDHeNpG63S7vdptVq0W636Xa7\nVZckrSpL+SiktKK63S4TExPMzMzsaOsvj42NVVWWtKo4oaqR0263mZ6eHto+NTVVQUXSSBs6oWrP\nXUNFDP37Uqler1d1CdKq4Zi7hsrMym7r1q0bWlOj0djLZ0FavQx3jZzx8XGazeYubc1mk06nU1FF\n0urjsIxGTn/SdHJykl6vR6PRoNPpOJkqLYITqpK0uvkNVUlaKwx3Saohw12Sashwl6QaMtwlqYYM\nd0mqIcNdkmrIcJekGjLcJamGDHdJqiHDXZJqyHCXpBoy3CWphgx3Saohw12Sashwl6QaMtwlqYYM\nd0mqIcNdkmrIcJekGjLcJamGDHdJqiHDXbvodru0221arRbtdptut1t1SZKWYN+qC9Do6Ha7TExM\nMDMzs6Otvzw2NlZVWZKWIDKz6hoAiIgclVrWqna7zfT09ND2qampCiqStAAxrNGe+wiKGPpnVZle\nr1d1CZIWyTH3EZSZldzWrVs3tJ5Go7GXz4CkPWW4a4fx8XGazeYubc1mk06nU1FFkpbKYRnt0J80\nnZycpNfr0Wg06HQ6TqZKq5ATqpK0ug2dpNujYZmI2BIR2wdumyPisIi4KSLujYhLB567PiKu2JP9\nSZIWZo967hHxDWAL8M6y6R7gF4CXAe8FLgKOAX4E3ASckJk/2M227LlL0uItf8+93OgW4PrMvCYz\nPwkcANwNbCqf82jgj4CNuwt2SdLyWo6e+5EUIf894PeB2yiCvUHRW38N8H7guMx8aI5t2XOXpMUb\n2nPf03A/H/gKsD/wNuAwimGYByhC/xbg48ClwFHAuRT/CZydmbfP2lZedNFFO9ZbrRatVmvJtUnS\nGrH84b7LhorJ0/VAuxyeISJOA/4b8GsU4+4t4BygkZlnznq9PXdJWrzl/fmBiDgOeCtFz3xf4Czg\nfuDW8vH9gIuB04F9ygJeBjwTuGOp+5UkzW9PvsT0XYoJ2T+kmES9DXhzZt5dPj4BfCIzvwoQERcA\n51FMtl6yB/uVJM3DLzFJ0uq2Ih+FlCSNIMNdkmrIcJekGjLcJamGDHdJqiHDXZJqyHCXpBoy3CWp\nhgx3Saohw12Sashwl6QaMtwlqYYM9xHS7XZpt9u0Wi3a7TbdbrfqkiStUnvyk79aRt1ul4mJCWZm\nZna09ZfHxsaqKkvSKuVP/o6IdrvN9PT00PapqakKKpK0SizvlZjqKGLoOapUr9erugRJq5Bj7gMy\ns7LbunXrhtbUaDT28lmQVAeG+4gYHx+n2Wzu0tZsNul0OhVVJGk1c1hmRPQnTScnJ+n1ejQaDTqd\njpOpkpbECVVJWt28hqokrRWGuyTVkOEuSTVkuEtSDRnuklRDhrsk1ZDhLkk1ZLhLUg0Z7pJUQ4a7\nJNWQ4S5JNWS44+XtJNXPmv9VSC9vJ6mO1vyvQnp5O0mr3OhfZm+ULnPn5e0krWYjNebu5e0kaXmM\nVLhXwcvbSaqjkRqWqYKXt5NUR2t+QlWSVrnqLrMXEc+LiC9GRC8ivhARx++N/UrSWrXi4R4RDeDD\nwCOB1wOPBz4UEWt+vH93brjhhqpLGBmei508Fzt5Lua3NwL2RcAhwDsy83LgKuAooLUX9r0q+Rd3\nJ8/FTp6LnTwX89sb4X5UeX/XrPujhjxXkrQMqhgaGZ1vKklSTa34p2Ui4jTgI8DvZebbI+ItwAXA\nKZm5aeB5flRGkpYgM3+q07w3wn1/4JvA/cDbKYK9BzzFzz5K0spY8WGZzNwGnAHcB/wJcDdwhsEu\nSStnZL7EJElaPpV/1twvOBUi4piI2BQR34+IeyNiOiKOrrquKkVEIyK+GhHbI2Ky6nqqEhEHRsR7\nI+KHEfHjiLix6pr2toh4cvn3YHtEXDDQflW/vcr6RlGl4e4XnHZxWHn/B8C7gRcCf1ldOSPhD4DD\ny+W1/BbzXcCZwJXABHBHteVU7jcBIuJRwEvLtrX892Ooqn84rP8Fpzdm5uURcShwIcUXnD5dZWEV\n+Gxm/kJ/JSJeATytwnoqFRHHUfyHfyHFRPyaVL57Ow14P3A+sD0z31VtVZX6OnB0RLSAoyky7C52\ndo5UqrqH7BecSpn57/3liDgBeCzwmeoqqk75zu0vgT8DPl9xOVXr/wf/HODfgPsi4m0V1lO124HP\nAa8CzgauBX5YaUUjqupwn23Nf8EpIo4F/hr4BrBWf1T+bOBJwPuAI8q2AyPi4OpKqsz+5f0BFEMQ\n/wC8KSJOqa6kSiXFMNVLgeeVy2s+N4apOty/Xt4/sbw/fFb7mhIRTwNuoPgewC9m5tZqK6rMEcDj\ngFsoAh7gFcDFlVVUnf6/hb/LzGuBD5bra3my/QPAg8C3gb+tuJaRVfWY+8eB7wK/ExH3AedQ9Fhv\nqLKoKkTEE4FNwEHAZcCJEXFiZn6g2soqcQ1wa7n8dGADxd+Vd1ZVUFUyc3NE3Aq8MCJeTfGu5kGK\nHvyalJk/johXAfdmZpbXXrb3Pkul4Z6Z2yLiDODPKb7g9CXg1Wv0C05Nit5qApeUbUnRS1lTMvN2\nirFVIuIHZfNMZm6urqpK/TrFHMRGim97n5WZX662pGpl5jWDq/hpmZ/il5gkqYaqHnOXJK0Aw12S\nashwl6QaMtwlqYYMd0mqIcNdkmrIcJekGjLcJamG/j9XKHW3dU1MsQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xc49b438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(0,10,1)\n",
    "y1 = np.array([1,1,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan])\n",
    "y2 =  np.array([np.nan,3,3,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan])\n",
    "y3 =  np.array([np.nan,np.nan,6,6,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan])\n",
    "y4 =  np.array([np.nan,np.nan,np.nan,9,9,np.nan,np.nan,np.nan,np.nan,np.nan])\n",
    "y5 =  np.array([np.nan,np.nan,np.nan,np.nan,12,12,np.nan,np.nan,np.nan,np.nan])\n",
    "y6 =  np.array([np.nan,np.nan,np.nan,np.nan,np.nan,15,15,np.nan,np.nan,np.nan])\n",
    "y7 =  np.array([np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,18,18,np.nan,np.nan])\n",
    "y8 =  np.array([np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,np.nan,21,21,np.nan])\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x,y1,'k')\n",
    "ax.plot(x,y2,'k')\n",
    "ax.plot(x,y3,'k')\n",
    "ax.plot(x,y4,'k')\n",
    "ax.plot(x,y5,'k')\n",
    "ax.plot(x,y6,'k')\n",
    "ax.plot(x,y7,'k')\n",
    "ax.plot(x,y8,'k')\n",
    "plt.xlim(0,10)\n",
    "plt.ylim(0,23)\n",
    "ax.plot([1,2,3,4,5,6,7,8],[1,3,6,9,12,15,18,21], 'ko',markerfacecoloralt='w')\n",
    "ax.spines['top'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.set_title('Average Service Rate:$\\mu(l)$')\n",
    "plt.xticks((0,2,4,6,8),('0','2','4','6','M'))\n",
    "plt.yticks((0,5,10,15,20), ('0','5%','10%','15%','$\\mu(k)$'))\n",
    "plt.savefig(r'D:\\average_service_rate')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 545,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "l"
      ]
     },
     "execution_count": 545,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import *\n",
    "l, p_0 = symbols('l, p_0_0')\n",
    "var('l')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 546,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.99468085,  0.85      ,  0.748     ,  0.68      ])"
      ]
     },
     "execution_count": 546,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "M = 4\n",
    "k = np.arange(0,M+1,1)\n",
    "k_l = k*l\n",
    "mu_k = np.where(k<2,0.188,np.where(k<3,0.220,np.where(k<4,0.250,np.where(k<5,0.275,np.where(k<6,0.300,np.nan)))))\n",
    "rou_k = 0.187/mu_k[1:]\n",
    "rou_k  # k >= 1\n",
    "mu_k\n",
    "rou_k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 547,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#from sympy import *\n",
    "#var('n,x')\n",
    "#Sum(n**2,(n,1,2)).evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 548,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 548,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from functools import reduce\n",
    "reduce(lambda a,b: a*b, range(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 549,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 549,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "l_rou_k = (l*(1-rou_k[:-1]**l)*rou_k[:-1])/(1-rou_k[:-1])*p_0\n",
    "l_rou_k.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 550,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.34388843085106*l*p_0_0], dtype=object)"
      ]
     },
     "execution_count": 550,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "k_equal_M = np.array([l*rou_k[-1]/(1-rou_k[-1])*reduce(lambda a,b: a*b, rou_k[:-1])*p_0])\n",
    "k_equal_M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 551,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99468085106382975, 0.845478723404255**l,\n",
       "       0.748**l*(0.845478723404255**l)**l], dtype=object)"
      ]
     },
     "execution_count": 551,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_rou_pow():\n",
    "    lst =[]\n",
    "    for i in range(1,M+1):\n",
    "        lst.append(reduce(lambda a,b: a**l*b**l, rou_k[:-1][:i]))\n",
    "    return np.array(lst[:-1])\n",
    "get_rou_pow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 552,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4,)"
      ]
     },
     "execution_count": 552,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_k_dot = np.concatenate((get_rou_pow()*l_rou_k, k_equal_M))\n",
    "p_k_dot.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 553,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.742063492063492*0.748**l*l*p_0_0*(-0.748**l + 1)*(0.845478723404255**l)**l + 1.24666666666667*0.845478723404255**l*l*p_0_0*(-0.85**l + 1) + 34.9689999999998*l*p_0_0*(-0.99468085106383**l + 1) + 0.369569318484042*l*p_0_0"
      ]
     },
     "execution_count": 553,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_L = reduce(lambda a,b: a+b, mu_k[1:]*p_k_dot)\n",
    "mu_L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 554,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.277531746031746*0.748**l*l*p_0_0*(-0.748**l + 1)*(0.845478723404255**l)**l + 0.466253333333333*0.845478723404255**l*l*p_0_0*(-0.85**l + 1) + 13.0784059999999*l*p_0_0*(-0.99468085106383**l + 1) + 0.138218925113032*l*p_0_0 - 0.034969)/((0.742063492063492*0.748**l*l*p_0_0*(-0.748**l + 1)*(0.845478723404255**l)**l + 1.24666666666667*0.845478723404255**l*l*p_0_0*(-0.85**l + 1) + 34.9689999999998*l*p_0_0*(-0.99468085106383**l + 1) + 0.369569318484042*l*p_0_0)**2*(0.742063492063492*0.748**l*l*p_0_0*(-0.748**l + 1)*(0.845478723404255**l)**l + 1.24666666666667*0.845478723404255**l*l*p_0_0*(-0.85**l + 1) + 34.9689999999998*l*p_0_0*(-0.99468085106383**l + 1) + 0.369569318484042*l*p_0_0 - 0.187)**2)"
      ]
     },
     "execution_count": 554,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_L = (0.187*(2*mu_L-0.187))/(mu_L**2*(mu_L-0.187)**2)\n",
    "D_L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 555,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_k_mole = get_rou_pow()*0.187"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 556,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99468085106382975, 0.845478723404255**l,\n",
       "       0.748**l*(0.845478723404255**l)**l,\n",
       "       0.68**l*(0.748**l*(0.845478723404255**l)**l)**l], dtype=object)"
      ]
     },
     "execution_count": 556,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_rou_pow_M():\n",
    "    lst =[]\n",
    "    for i in range(1,M+1):\n",
    "        lst.append(reduce(lambda a,b: a**l*b**l, rou_k[:i]))\n",
    "    return np.array(lst)\n",
    "get_rou_pow_M()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 557,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.125*0.748**l*l*(0.845478723404255**l)**l + 3.25527297419647*l,\n",
       "       2.125*0.748**l*l*(0.845478723404255**l)**l + 3.27268085106383*0.845478723404255**l*l,\n",
       "       5.39768085106383*0.748**l*l*(0.845478723404255**l)**l,\n",
       "       3.27268085106383*0.68**l*l*(0.748**l*(0.845478723404255**l)**l)**l + 2.125*0.748**l*l*(0.845478723404255**l)**l], dtype=object)"
      ]
     },
     "execution_count": 557,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nol_rou_k = ((1-rou_k)*rou_k)/(1-rou_k)\n",
    "r_k_deno = nol_rou_k.sum()*l*get_rou_pow_M()+l*rou_k[-1]/(1-rou_k[-1])*reduce(lambda a,b: a**l*b**l, rou_k[:-1])\n",
    "r_k_deno"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 558,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (3,) (4,) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-558-66fa10fceb38>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mr_k\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mr_k_mole\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mr_k_deno\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mr_k\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (3,) (4,) "
     ]
    }
   ],
   "source": [
    "r_k = r_k_mole/r_k_deno\n",
    "r_k.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 559,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.656910908511448*0.748**l*l*p_0_0*(-0.748**l + 1)*(0.845478723404255**l)**l + 1.57201794770945*0.845478723404255**l*l*p_0_0*(-0.85**l + 1) + 67.8030888471351*l*p_0_0*(-0.99468085106383**l + 1) + 0.250792417040337*l*p_0_0)/l**2"
      ]
     },
     "execution_count": 559,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "average_N_t = reduce(lambda a,b: a+b, p_k_dot*(rou_k+np.e**(-rou_k)-1))/l**2\n",
    "average_N_t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 560,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a, b, e, h = symbols('a, b, e, h')\n",
    "F_l = (a*mu_L+ b*average_N_t + 100*e*r_k.mean()+ h*D_L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 561,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "F_L = F_l\n",
    "lst_result = []\n",
    "for i in range(1,101):\n",
    "    result = F_L.subs([(p_0, 0.01), (l, i), (a, np.e), (b, 10), (e, np.e), (h, 100)])\n",
    "    lst_result.append(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 568,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApAAAAGJCAYAAAA9hAOeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8XGV97/HPb+dGIISbJ0guCEhiRcSiqVDBsr3UHrBW\nLE1LvQBKi9WjgkqFUyhEFBRRi9RSrTeqtVjk4pWq9CBSEUURBcFyDYiBECBASEKSfXnOH2tNMhlm\nz57Ze6+5rPm8X6/1WjPrWWvNs2fB3t88v3WJlBKSJElSswY63QFJkiT1FgOkJEmSWmKAlCRJUksM\nkJIkSWqJAVKSJEktMUBKkiSpJQZIST0tIi6KiNGIOHMK9zk93+/j+b7Pm6p9T6AvU/7zTbAfgxFx\nW0QM5/3ZfoL7WRIRQxHxxfz9Xvn+RqvW+WFErImInaeq/5KmlgFSUteKiHsr4WKM6bCq1afyprZ/\nChwDbAYuAH4whfuuq16Qyn0XOB+4vug+jONC4HfY2p+h2hUi4poxjtNrq1Z7P9nfno/XbF59/D4O\n7AycPKU/gaQpM73THZCkBj4H7JK/fjswE7gU+G2+7Lf1NpoCS/L5lSmlkwr6jLFsE4RTShcDF7e5\nD/UsIevbO1JKK8ZZ9wfATVXv7wKIiHnAUcAdKaVfNNj+W8A64PiIODOlNDLxbksqgiOQkrpWSukD\nKaX3pJTeA2wkCzCfrCxLKd1dtfozIuLyiNgQEb+MiBdUGiJiz4j4SkSsjIjHIuK7EfG8ep8ZEcuB\ns/K3x+QjaMfWlpLHKL1WRtzeERF3RMSTEfGliJhRtc4ReYn2sYh4IiK+FRHPAu7ZusqW/TyrXgk7\nIk6IiFsiYl1E3BkRH4iIWXnbYL7+ioj4u4hYHREPRUTD0bxx9jlK9vcigLsjYrwAeUXVMXpPSunW\nfPmryAYuGo7oppQ2Az8BdgeWjvNZkjrAACmpLP4PWcC8B3g+8I8A+bl6VwPLgF8A3wAGgasjYrc6\n+7meLLwA3EZWrr2tqr22VF6vdL4c+CHZ79g3AG/K+/KHZKNrLwGuIxtNXQSsBb5Qtf35+bS29nMi\n4u3Ap4AFwFfIAtlpwCdq+vAs4PXAtcD/As6NiH3r9LWZfVbv+/NkI8ON/GlEnF81VUaRD8jnvx5n\n++p1freJdSW1mQFSUllcmVI6CnhH/r4SPF4N7AM8ANwBPAbcTxaq/qx2Jyml75Kd5wdwQz6C9tMW\n+/LWlNJbgK/W9OXEfP6JlNIfp5SOB16YUnqMraOeqWrk7rE6+678fCemlP4KqJxfeHxlxDA3DLw8\npfRnZD9vsDXAtbLPmSmld1f6BpyVUvpgox8e+APgXfn0TmDHfHklSD45zvawNTzv0nAtSR3hOZCS\nyqJyzt0T+XyHfL5XPl/A1gCX8unZk/i8aWMsT3X6Mief753Pf7xl5a3n90WTn7tX/hmVEbrb8/kA\nsLBqvVUppdX568fJRjrnUF+jfS4C7q6zTSMnpZQuqLP88Xy+Y522WnNrtpHURRyBlFQWw/m8tqRc\nOV/vZymlgXyaBuwKnNPC/tfn853y+f4T6EvlPMeDKwsiovIP+ZGqZY3C5AqysPnc/P1z8vko2Uhj\nbR/q9WOi+5ysX+bz/ZpYt9KXmxquJakjDJCSekmzo3TVriQLSEsj4kcR8amIuJKspD1WSbeen+fz\nYyPiI8A/jdGnRn2sjMqdGBHfjIjPATfmyx4iu21QABdHxIfH2Eflcz+Rb//1/P3n8otPJqLZfU7k\n+6/2PbLb//xBo5XyUvxBZN9Jq6cPSGoDA6SkXlEpOze7PGtMaQPwcrJb4Swiu7/jYuBLbC3VNrPP\nLwH/TnbqzxFsvY9ho4tqttlPSukqsnMyfwQcSnZLm9/kbZuBU4CHgT8nu21RvX1cmLetBI4mG2k8\nh23L8838PFsbx99nvZ+t7q7G+ZzVwGXAcyLiwAb7+WOycvtnU0q198WU1AUipam8926LHx5xAfAX\nZCezfzul9Jp8+SHAP5Pdd+xW4K9SSjd1W5skqTURsZjsd+lXUkrHjLHOf5OVufdJKT1Rbx1JndXp\nEcjE1hvkVm5RsR3Zv1B3AE4iuw/YpZHplrZOf2+S1JNSSnemlGaOFR7zdV6aUtrN8Ch1r45ehZ1S\nOjG/ge67qhYfDswD/jal9KmI2AP4e+BlZCevd0PbINl95SRJkvpON9zGp/ak7MptLlbWzPdh620d\nOt1W6aMkSVLf6YVSbCVg1jtZs1NtkiRJfasbRiBrVe7ZtiifL8jn97D1/mvd0LaNiEhnnrnlUbUM\nDg4yODhYu5okSVI3mNSgWKevwn412c14PwTcTPbs2huAq4ANwHnA6cBGYF9gJnBfN7Slmi8uImoX\nSZIkdatJBchOl7BPJguPieyGvv8CHAgsA9YB5wOrgGUps6lb2gr+XiRJkrpWR0cgy8QRSEmS1EN6\negRSkiRJPcYAKUmSpJYYICVJktQSA6QkSZJaYoCUJElSSwyQkiRJaokBUpIkSS0xQEqSJKklBkhJ\nkiS1xAApSZKklhggJUmS1BIDpCRJklpigJQkSVJLDJCSJElqiQFSkiRJLTFASpIkqSUGSEmSJLXE\nAClJkqSWGCAlSZLUEgOkJEmSWmKAlCRJUksMkJIkSWqJAVKSJEktMUBKkiSpJQZISZIktcQAKUmS\npJYYICVJktQSA6QkSZJaYoCUJElSSwyQkiRJaokBUpIkSS0xQEqSJKklBkhJkiS1xAApSZKklhgg\nJUmS1BIDpCRJklpigJQkSVJLDJCSJElqiQFSkiRJLTFASpIkqSUGSEmSJLXEAClJkqSWGCAlSZLU\nEgOkJEmSWmKAlCRJUksMkJIkSWqJAVKSJEktMUBqi6uvvpoLL7yw092QJEldLlJKne5DKURE6vXv\nct999+Xuu+/mnnvuYe+99+50dyRJUnFiMhs7AikARkdHue+++wBYtWpVh3sjSZK6mQFSADz88MMM\nDw8DsGbNmg73RpIkdTMDpAB44IEHtrx+9NFHO9gTSZLU7QyQAmDlypVbXhsgJUlSIwZIAdsGSEvY\nkiSpEQOkAEcgJUlS8wyQAjwHUpIkNc8AKcAStiRJap4BUoAlbEmS1LyuDZARcW9EjFZNN+XLD4mI\nmyNiY0TcGBEHVm3T1rYyMUBKkqRmde2jDCNiBXAv8M/5oseA/86XrQfOA04HNgH7ArPa2LY4pTRa\n09+efZThxo0bmT179pb3c+bM4cknn+xgjyRJUsEm9SjD6VPViwIEWXi7MqW0DiAiXgfMA/42pfSp\niNgD+HvgZcBObWwbBK5uw3fQFpULaBYtWsQDDzzAunXr2Lx5MzNnzuxwzyRJUjfq2hI2kIBjgLUR\n8VBEvAXYO29bWTPfp81tleWlUClfL1y4kF133RWwjC1JksbWzQHyM8Ay4A3ARuDTddapDL/Wqx0X\n2VYqlQA5f/58dtttN8ArsSVJ0ti6toSdUjqn8joilgLvBn6bL1qUzxfk83vIys3tbHua5cuXb3k9\nODjI4OBgvdW6TqWEvWDBAlatWgU4AilJksbWlQEyIg4Azgb+k6yPxwAbyC6iWQ28LSLWAccDK4Br\ngJltbnua6gDZSyojkAsWLGDFihWAAVKSJI2tW0vYq8n69n7gQ2Sh7XUppQfJytrrgPOBVcCylNnU\nzrZ2fAntUh0gLWFLkqTxdOUIZEppFfDqMdr+GzigG9rKot45kI5ASpKksXTrCKTaqPocSK/CliRJ\n4zFA9rmUkiOQkiSpJQbIPrdmzRo2bdrE3LlzmTNnjudASpKkcRkg+1z1BTSAI5CSJGlcBsg+V33+\nI+A5kJIkaVwGyD431gikJWxJkjQWA2Sfq76ABrYtYZfsdpeSJGmKGCD7XO0I5OzZs9luu+3YvHkz\n69ev72TXJElSlzJA9rnacyDBMrYkSWrMANnnakcgwSuxJUlSYwbIPld7DiR4JbYkSWrMANnHNm/e\nzOrVqxkYGGD33XffstwStiRJasQA2cdWrVoFwDOf+UymT5++ZbklbEmS1IgBso/VO/8RDJCSJKkx\nA2Qfq3f+I2w9B9IStiRJqscA2cccgZQkSRNhgOxj9e4BCQZISZLUmAGyj401AmkJW5IkNWKA7GNj\nnQPpCKQkSWrEANnHPAdSkiRNRKSUOt2HUoiI1EvfZUqJOXPmsGHDBh5//HF22mmnLW1DQ0PMnDmT\ngYEBhoaGGBjw3xmSJJVMTGZjk0GfWrt2LRs2bGCHHXZg7ty527TNmDGDuXPnMjo6yhNPPNGhHkqS\npG5lgOxT1ec/Rjz9HyGWsSVJ0lgMkH1qrPMfKypXYhsgJUlSLQNknxovQFZGIL2VjyRJqmWA7FOV\nm4jX3sKnwhK2JEkaiwGyTzU7AmmAlCRJtQyQfarZcyAtYUuSpFoGyD7lCKQkSZooA2Sf8hxISZI0\nUQbIPjQ8PMyqVasA2GOPPequYwlbkiSNxQDZhx566CFGR0eZN28eM2fOrLuOI5CSJGksBsg+NN75\nj2CAlCRJYzNA9qHxzn8ES9iSJGlsBsg+dP/99wONRyB32mknBgYGWLt2LUNDQ+3qmiRJ6gEGyD70\n05/+FID9999/zHUGBgYchZQkSXUZIPvQddddB8Chhx7acD2fhy1JkuoxQPaZBx98kHvuuYc5c+bw\n/Oc/v+G6lRFIL6SRJEnVDJB9pjL6ePDBBzN9+vSG63oltiRJqscA2WeaLV+DJWxJklSfAbLPVALk\nIYccMu66lrAlSVI9Bsg+sn79em666SYGBgY46KCDxl3fErYkSarHANlHbrjhBoaHh3nBC17Ajjvu\nOO76BkhJklSPAbKPtFK+Bp9GI0mS6jNA9pFWLqABRyAlSVJ9Bsg+MTIywvXXXw80PwJpgJQkSfUY\nIPvErbfeyhNPPMGee+7JwoULm9rGErYkSarHANknWj3/ERyBlCRJ9TV+FAkQEfOBlwB75YvuA36Y\nUnqwwH5pirV6/iPA9ttvz6xZs9i4cSMbNmxg++23L6p7kiSph4wZICPitcB7gEOBqGlOEXEt8A8p\npW8U2D9NkYmMQEYEu+22Gw888ABr1qwxQEqSJKBxCfuKvP3vgFcC+wHPA/4wXzYD+FrRHdTkrVy5\nknvvvZe5c+ey//77t7StT6ORJEm1GpWw908p3VZn+a+B/wecGxH7FdMtTaXK6OPBBx/MtGnTWtrW\n8yAlSVKtRgFyVUTsOlZjSmnNGAFTXWYi5z9WVAKkV2JLkqSKRgHyESDVWR758taGstQxEzn/scIS\ntiRJqtUoQF7boK1esFQXWrduHb/4xS+YNm0aBx10UMvbW8KWJEm1xgyQKaXBNvZDBfnJT37CyMgI\nL3rRi9hhhx1a3t4StiRJqjXmVdgR8daIGPO+LRGxQ0S8tZhuaapMpnwNlrAlSdLTNSph/wPw8Yj4\nPvBT4AGy8x/nA78HvIyslP3pojupiZvMBTSwdQTymmuu4ec//zkvfOELp6xvkiSpNzW6D+TewMeB\n5wJnkgXFTwFnAM8BzsvXKb2IOCQibo6IjRFxY0Qc2Ok+jWdkZITPfvazXHttdirrREcgX/GKV7B4\n8WLuvfdeXvziF3PKKafw1FNPTWVXJUlSj4mUxr8eJiL2APbM3/6mnx5jGBHbAfcC68lC8+nAJmBx\nSmm0ar3UzHfZDtdeey0nnXQSN910EwDLli3jkksumfD+1q9fzxlnnMH555/P6Ogoixcv5jOf+QyH\nHXbYVHVZkiS1V+1TBlvbeKzQExFzgQ0ppeHJfECvi4jXAZcBf5tS+lhEvB/4e+CVKaWrq9breIBc\nsWIF73vf+7j00ksBWLhwIR/5yEc4+uijiZjUfydAdkHO8ccfz6233grAcccdxyte8QqWLFnCkiVL\n2HnnnSf9GeoPJ77+RH721Z8xK81iU2xi6bKlfOLfP9HpbklSPyksQI4CRwNXAT8HXp9Sun4yH9aL\nIuI9wEfJfv6vRMQJZKX8v04pfa5qvY4EyCeffJLvfOc7XHHFFVx++eVs2rSJ2bNnc+qpp3LyySdP\n+fOrN2/ezIc+9CHOPvtshoaGtmmbN28eS5YsYf78+cydO3ebaccdd2TWrFnMnDlzm2nGjBlMmzZt\nm2n69OkMDAwQEQwMDGx5XZlq31fCcb33FfWWVau3bqP1x9vHRExFyO8Fp7/1dFZesZKTOXnLso/y\nURa8bgEf/PQHO9gzSe3W6YGXfjZv3rxJ/dFpdBFNxTTgWcDsyXxQiXT8r/zDDz/M17/+db72ta/x\nX//1X2zatGlL2xve8AY+/OEPs3DhwkI+e+bMmZx55pksW7aMiy++mNtvv5077riDO+64g9WrV7N6\n9epCPlflsZSlnMd52yw7mZN53xXvY94V8zrUK0nqL5MN780EyH7/58E9+XxRPl9Qs3yL5cuXb3k9\nODjI4ODglHUipcSPf/xjLrjgAi699FKGh7MzCyKCQw45hCOPPJIjjzySfffdd8o+s5H99tuPD3zg\nA1vej46OsnLlSm6//XYeeeQR1q5du8305JNPsnnz5rrTyMjINtPw8DApJUZHR7fMK68rU/X7yvdT\n+776u6tdVq3euo3WH28fE9FP/wrf/pH6o+Kzmc0znvGMNvdGUqf1S/WlbMYrYd8PrCO7EnsF2YUk\nAKSUDmhHBzstImYB9wEb2HoRzUZg3+qadVEl7E2bNvEf//EfXHDBBdx4440ATJs2jVe96lW87nWv\n40/+5E/Yfffdp/xzpaIcMuMQzh4++2nLT5t+GtcNXdeBHklSXyq0hL2o6nVf3LKnVkppU0QsA/4J\nOB/4Fdn5j4UPGX3+85/n1FNP5eGHHwayezKecMIJvO1tb2PRokXjbC11p6XLlvLRiz+6zTmQ53Ee\nS5ct7WCvJEmtaOo2PhrfVI9Ann/++bz73e8G4Hd/93d517vexdFHH83s2Z6Kqt7nVdiS1HHFXIWt\n1kxlgPzYxz7GySdnozOf/OQnefvb3+45IpIkaSoVfhW22ugjH/kIp5xyCgCf/vSnOeGEEzrcI0mS\npG01epSh2uycc87hlFNOISL43Oc+Z3iUJEldacwAGRFfjIiDIuKMiHheOzvVj8466yxOO+00IoIv\nfOELvOUtb+l0lyRJkupqdBufEeBtZE9dOTqlNPGHKfeByZwD+b3vfY8/+qM/YmBggIsuuog3velN\nU9w7SZKkbUzqHMhGJewHycIjwFciYiSfRvNwqSly2WWXAXDqqacaHiVJUtdrdBHNe4GTgRcBvwXW\nVrV56fYUSSlx5ZVXAnDUUUd1uDeSJEnjG/c2PhFxEXBhSumGtvSoR020hH3LLbdwwAEH8MxnPpOV\nK1cyMOB1TZIkqXCFlbABSCkdBzwvIv4jn46dzAdqW5XRx8MPP9zwKEmSesK494GMiNOBs6oWLYuI\nhSmlpz/MVi2rBMgjjjiiwz2RJElqTjMl7HuBm4H35Is+BrwgpbRXoT3rMRMpYT/++OM84xnPAODR\nRx9lp512KqJrkiRJtYotYQO7AN9LKd2VUroL+C9g18l8qDJXXXUVIyMjHHrooYZHSZLUM5p5lOHP\ngHMi4qD8/WuBnxbXpf5h+VqSJPWiZkrYzwO+AeydL7obeG1K6baC+9ZTWi1hj46OMn/+fB566CFu\nueUW9t9//wJ7J0mStI1JlbDHDZAAETEDeE7+9n9SSsOT+dAyajVA3njjjSxdupRFixZx3333ETGp\n4yhJktSKSQWPZkrYpJSGgF9N5oO0reryteFRkiT1Em882CGe/yhJknpVUyVsja+VEvYjjzzCvHnz\nmDFjBo8++ihz5swpuHeSJEnbKPY2PhHxjIjYPX/9ioh4Y0RsN5kP7Xff/e53SSlx2GGHGR4lSVLP\naeYcyG8BN0XEJcBV+bIjgNcX1quSs3wtSZJ6WTO38XkCeC/wbOBQ4DZgWUrJm4lXabaEPTIywrx5\n81izZg233347S5YsaUPvJEmStlH4k2gGgAVk4fE7wHWAJewJuuGGG1izZg3PfvazWbx4cae7I0mS\n1LJmAuQNwJnAIWQl7MXAiiI7VWbevkeSJPW6Zs6B/EvgDcAdKaUbImJP4EfFdqu8PP9RkiT1unFH\nIFNKq4GrgcUR8WxgNXBr0R0rq7vuuguA3/u93+twTyRJkiZm3BHIiDga+Deyky1vAU4F1gNHFtu1\ncnrqqacA2GGHHTrcE0mSpIlp5hzI95ONQFZO2Ps28JLCelRiIyMjDA0NERHMmjWr092RJEmakGYC\n5HyyAAmQgCFgdmE9KrGNGzcCsN1223kBjSRJ6lnNXETzK+CY/PWbgMOBXxbWoxKrlK+32867IEmS\npN7VzAjke4Dd89fHkoXOkwvrUYlVRiBnz3YAV5Ik9a5xRyBTStdHxGLg98lK2NenlB4rvGclVBmB\nNEBKkqRe1sxV2MeSBceK10QEKaUvFtetcrKELUmSyqCZcyC/UGdZAgyQLbKELUmSyqCZAPm+qtc7\nk50H+cNiulNulrAlSVIZNHMO5Eer30fEzcDphfWoxCxhS5KkMmjmHMhvsvUcyOnAUmBGkZ0qK0vY\nkiSpDJopYb+65v1GsscZqkWWsCVJUhk0EyD3qXo9DDyUUhoqqD+lZglbkiSVwZgBMiKOYtvb91S3\nkVK6vLBelZQlbEmSVAaNRiC/2qAtAdOmuC+lZwlbkiSVQaMAeVbbetEnKiOQlrAlSVIvGzNAppSW\nt7EffcERSEmSVAbN3MZnPrAc2B/YMnSWUnphcd0qJy+ikSRJZdDMVdifAf4IGACGyO4B+USRnSor\nL6KRJEllMNDEOi8BPpy/fjXwaeCfCutRiVnCliRJZdBMgJwJrMhfHwisBd5VWI9KzBK2JEkqg2ZK\n2PcBuwG3AOfmy/6nsB6VmCVsSZJUBs0EyGXAZuA/gdPI7gF5dpGdKitL2JIkqQyaCZB/DfxbSuln\nwNEF96fULGFLkqQyaOYcyHcBN0TEryPi9IjYq9gulZclbEmSVAbNBMj9gfcDw2RPp7k7In5YaK9K\nyhK2JEkqg3EDZErptpTS+8nuBVm5fc9LCu1VSVnCliRJZdDMk2jeDfwZcDAQwN3Avxfcr1KyhC1J\nksqgmYtoPgY8DFxIdjHNT4rtUnlZwpYkSWXQTID8Y+C7KaWRojtTdpawJUlSGURKqdN9KIWISI2+\ny5QS06dPZ3R0lKGhIaZPbya7S5IkFSIms3EzV2FrCgwNDTE6Osr06dMNj5IkqacZINvE8rUkSSqL\npgJkRDwzIo6MiPkRsWdEzC2qQxGxPCJGa6YD8rYdI+LiiFgfEQ9GxHurtmtrW6u8AluSJJXFuAEy\nIl4J3AlcBjwXuBT4VMH9guyxiZXpvnzZB4G/AM4FrgfOi4iXdaitJV6BLUmSymLci2gi4hfAZmAp\n8IfAC4GTUkoLCulQxHLgDGAHYHP11d8R8Thwf0rp+RGxN9k9Kb+UUjq2jW3/llI6pk6/G15E8+tf\n/5r99tuPJUuWcPvtt0/FVyVJkjRRhV9Esy/Z6CNAAh4Ddp7MhzZpHbAhIr4SEbMjYldgLrAyb6/M\n94mIXdrYtvdEfhhL2JIkqSyauRz4buDI/PUfAkcBkxpCi4jfAvPrNL0ZuBE4AXggf//nwK/Y+hjF\nLbtp9BFtbhuXJWxJklQWzQTI09g6AnkKWTn7Tyf5uS8FZtRZviqltLbyJiLuJwusz00pPRYRTwIL\n8+ZKCf2edreN9UMtX758y+vBwUEGBwe3vPcqbEmSVBbjBsiU0rci4vlko48JuCqldOdkPjSltGKs\ntoj4KnAzcD/wxnxx5fGJFwHvjIgzgAPz/lzUobanqQ6QtSxhS5Kksmj2jtbzgCfIzpn8/Yj4/ZTS\nFwvq023AcWQl7oeAc4B/zNtOB3YnGwldC5yaUvp+h9paYglbkiSVxbgBMiK+DPxlzeIEFBIgU0pn\nAmeO0fYk2W19Ot7WKkvYkiSpLJoZgXwN8DPgcmC42O6UlyVsSZJUFs0EyGuA61JK5xbcl1KzhC1J\nkspizAAZEd8kK1XvDJwdEa8B1lTaU0p/Unz3ysMStiRJKotGI5Cvrnn/kiI7UnaWsCVJUlk0CpCV\nJ65M6gbayljCliRJZdEoQF4DvCOl9O029aXULGFLkqSyaPQs7GcBO7SrI2VnCVuSJJXFeFdhHxYR\ndYfMCryReClZwpYkSWUxXoB8Wz7VKuxG4mVlCVuSJJXFeAHyy8Av6yxPBfSl1CxhS5KkshgvQH4z\npXRJW3pScpawJUlSWTS6iOY3wIZ2daTsLGFLkqSyGHMEMqW0Vxv7UXqWsCVJUlk0GoHUFLKELUmS\nysIA2SaWsCVJUlkYINvEErYkSSoLA2SbWMKWJEllYYBsE0vYkiSpLCIl7wk+FSIijfVdjo6OMm3a\ntC2vI6KdXZMkSao1qTDiCGQbVM5/3G677QyPkiSp5xkg28DytSRJKhMDZBt4BbYkSSoTA2QbeAW2\nJEkqEwNkG1jCliRJZWKAbANL2JIkqUwMkG1gCVuSJJWJAbINqm/jI0mS1OsMkG3gCKQkSSoTA2Qb\nGCAlSVKZGCDbwBK2JEkqEwNkGzgCKUmSysQA2QYGSEmSVCYGyDawhC1JksrEANkGjkBKkqQyMUC2\ngY8ylCRJZWKAbAMfZShJksrEANkGlrAlSVKZGCDbwBK2JEkqEwNkG1jCliRJZWKAbANL2JIkqUwM\nkG1gCVuSJJWJAbINLGFLkqQyMUC2gSVsSZJUJgbINrCELUmSysQA2QaWsCVJUpkYINvAErYkSSoT\nA2QbWMKWJEllYoBsA0vYkiSpTAyQBRsaGmJkZIRp06YxY8aMTndHkiRp0gyQBbN8LUmSysYAWTDL\n15IkqWwMkAXzCmxJklQ2BsiCWcKWJEllY4AsmCVsSZJUNgbIglnCliRJZWOALJglbEmSVDYGyIJZ\nwpYkSWVjgCyYJWxJklQ2HQmQEXF4RNwSEaP5tGtV244RcXFErI+IByPivd3Y1ixL2JIkqWw6NQI5\nG/gBcBeQato+CPwFcC5wPXBeRLysC9uaYglbkiSVTUcCZErp8pTSO4AH6jQfC9yaUjoLqIz4HddF\nbW9u9uflVSNmAAALsUlEQVQES9iSJKl8uuocyLyUPRdYmS+qzPeJiF26pG3vVn4mS9iSJKlsCguQ\nEfHbqnMcq6djWtlNj7SNyRK2JEkqm+kF7vulwIw6y1eNtUFKaU1ErAUW5osW5PN7UkqPRcST3dA2\nVv+XL1++5fXg4CCDg4OWsCVJUukUFiBTSivGaouIfYFBYA+ykb03RcSdKaUrgX8F3hkRZwAHkl1k\nc1G+6UVd1PY01QGywhK2JEkqmyJHIBs5FPgXskCWgH8ArgGuBE4HdgdOAdYCp6aUvp9v101tTbGE\nLUmSyqYjATKldBFjjOSllJ4Eju72tmZZwpYkSWXTVVdhl5ElbEmSVDYGyIJZwpYkSWVjgCyYJWxJ\nklQ2BsiCWcKWJEllY4AsmCVsSZJUNgbIglnCliRJZWOALJglbEmSVDYGyIJZwpYkSWVjgCyYJWxJ\nklQ2BsiCVUYgLWFLkqSyiJRSp/tQChGRar/LlBIDA1lGHxkZ2fJakiSpw2IyG5toClQZfZw1a5bh\nUZIklYappkCWryVJUhkZIAvkBTSSJKmMDJAFMkBKkqQyMkAWyBK2JEkqIwNkgRyBlCRJZWSALJAB\nUpIklZEBskCWsCVJUhkZIAvkCKQkSSojA2SBDJCSJKmMDJAFsoQtSZLKyABZIEcgJUlSGRkgC1QJ\nkI5ASpKkMjFAFqhSwnYEUpIklYkBskCWsCVJUhkZIAtkCVuSJJWRAbJAlrAlSVIZGSALZAlbkiSV\nkQGyQJawJUlSGRkgC2QJW5IklZEBskCWsCVJUhkZIAtkCVuSJJWRAbJAlrAlSVIZGSALZAlbkiSV\nkQGyQJawJUlSGRkgC2QJW5IklZEBskCWsCVJUhkZIAtkCVuSJJWRAbJAlrAlSVIZGSALMjw8zPDw\nMAMDA8yYMaPT3ZEkSZoyBsiCVJevI6LDvZEkSZo6BsiCWL6WJEllZYAsiFdgS5KksjJAFsQrsCVJ\nUlkZIAtiCVuSJJWVAbIglrAlSVJZGSALYglbkiSVlQGyIJawJUlSWRkgC2IJW5IklZUBsiCWsCVJ\nUlkZIAtiCVuSJJWVAbIglrAlSVJZGSALYglbkiSVlQGyIJawJUlSWRkgCzI8PExEGCAlSVLpREqp\n030ohYhItd9lSomUEgMD5nRJktRVYjIbT5+qXujpIoKISR0fSZKkrtORobGIODwibomI0Xzatapt\nedXyynRA3rZjRFwcEesj4sGIeG/Vdm1tkyRJ6ledqq3OBn4A3AWMVUM/umq6L1/2QeAvgHOB64Hz\nIuJlHWpTiVxzzTWd7oImyGPX2zx+vctj19siYnAy23ckQKaULk8pvQN4oMFq3wQuSyldklJ6Il92\nLHBrSuksoDIaeFyb297c5I+pHuIvwt7lsettHr/e5bHreYOT2bibr+5YB2yIiK9ExOy8zD0XWJm3\nV+b7RMQubWzbe0p+OkmSpB5VWICMiN/WOZdxNCKOGWfTG4ETgNcA3wD+nGz0r7bU3ejqlHa3SZIk\n9Y/KrWameiIbqVtSZ5pbtc41wAiw6xj7eD4wCnw5f/8E8Kv89T552792oq1OX5OTk5OTk5OTU69M\nk8l5hd3GJ6W0Yqy2iNiXrPa+B9nI3psi4s6U0pUR8VXgZuB+4I35Jj/J5xcB74yIM4ADyb6AizrU\nVvvzOkIpSZL6QkduJB4RxwGfJwtkkIXIa1JKL4+I95MFx/nAQ8C/AX+fUkoRsSPwGbLy9lrg4yml\n8/J9trVNkiSpX/kkGkmSJLWkm6/C7gkRcUhE3BwRGyPixog4sNN9Un0RsTgivh8Rj0TE2oj4XkTs\nk7d5HHtERGwXEbfnF+X9Y77M49flImLniPhiRDweEU9GxA/y5R67HhARp0bEb/LjdE9EvCNf7vHr\nMhFxQUQ8lP+O/GbV8jGP1USOowFyEiJiO+AyYAfgJGB34NKI8HvtTvPz+RnAF4BXAp+NiFl4HHvJ\nGcCC/HXy+PWMzwOvJzst6ETgTo9db4iI5wLnAJuBdwMzgAsiYiEev26UgIurXjfKKzHRLONBnpzD\ngXnAhSmlTwGfI7v6fLCTndKYfpRSellK6cKU0onAY8B+eBx7RmSPNT0JOLNq8RF4/LpaPtJ/JPDv\nwN+R3c3ir/DY9Yr1wDDZ/ZCvJrs+YSNwMB6/rpP/ffuHmsVj/Z17WYO2wUafU9hV2H2iclNxbzbe\nA1JKQ5XXEbEU2AW4FI9jT8j/NfxZ4JPAz6qa9srnHr/utV8+fzFZGBmJiE+QBRHw2HW1lNJvIuJv\ngH8Bfk12S7s3A3vmq3j8uk/tnWHG+ju3D9lDU+q1NTyOjkBOLW/l0wMi4nfIblK/AngnTz9uHsfu\n9GbgWcCXgIX5sp2BmTXrefy6z6x8vj3ZwyGuA97H0wcxPHZdKCL2BD4B/AJ4LfBLsn/I7VC7apu7\npomrHKt6V1I3dRwdgZyce/L5ony+oGa5ukxE7EdWgtkAvDyl9FBEeBx7w0Lgf5H98ap4I/5/2Asq\nx+K/U0pfi4jdgZez9Q+Vx667HUwWFi9PKX0zIl4AnEU2Ggkev15QuTd3vWO1U4O2MRkgJ+c/gdXA\n2yJiHXA82UG6ppOdUn0RsQj4PrAr2fkhvx8RBwNfw+PYCy4Bbslf7w8sJ/t/8Gzgcjx+XSuldFNE\n3AK8MiL+mmw0eRj4NvAePHbd7n/y+Zsi4iHgDWQjV3fg786uExGvJvsdCbBnRBwP3MDYx2pmg7ax\nFfUow36ZgJeSPTlnE9lzvF/Y6T45jXmsBsnO3RnJ56PAiMex9ybgsPz4XeDx642J7DzIHwFPkQWS\noz12vTMBbwfuyo/fXcDbPH7dOZENlFT/rRsBjml0rCZyHL2RuCRJklriRTSSJElqiQFSkiRJLTFA\nSpIkqSUGSEmSJLXEAClJkqSWGCAlSZLUEgOkJOUi4t6IGK0zrYiIZ+Wvv9npfgJExGDen3+cwn0+\nLyJuiYihiPhNnfZr8s/cdao+U1Jv8kk0krTVO8ie1/wasqdt/DPwA2B91TrddvPcqezPccDzgAuA\n/9eGz5PUoxyBlKRcSulbKaVL2Pq87Z+klC5JKX2brc9t3jkivhERayPiy5VtI+KPI+KXEbEuIn4R\nEa+s9xkRsTwiHoqIpyLizoj4y3z5iyPi2ny/D0TE6/LlP46IJyJifUT8LCIOHWO/z42Iq/J1742I\nk8ZYb5eIuCgiHo6I1RHxrxGxc0QcB7w3X+1dwLtb+/Yk9RNHICWpNS8BTgOeCfxlRPwz2XNkLwN+\nDnwAOAq4IiIWp5RWVTaMiF2AM4CrgX8F9skWx67AlWS/k5eTPXpsJN/se8BnyJ7hfiLweWBJdYci\nYhrwdbLR048ALwA+HhF3pZS+VdP/TwBvBM7L3/9tPj8z/6xXAWfh84wlNWCAlKTW/DildG5EBLAU\n2IsssM0ADsonyEq9BwNfq9p2HbAK+B3gUOAG4Arg5WQB8byU0scrK0fEHOBFwP8FplX2GxGzavr0\nHGDf/PUHqj7/lUBtgDwC+G1K6ZT8M94A/O+U0rERcRdZgLw6pXRts1+IpP5jgJSk1qzJ58P5fICt\n5wWeC1zF1nL3r6s3TCkNRcQLyEYoDwQ+BQwCF+erBNt6I3A48GXgi8CH8u1qA2TFd4CPVu1n1Rjr\nSdKkGCAlafKuAjaTBcO7gV2AZcCfVa+Ujyh+FLgOuBF4PbAH8CPgUeCtEfEQWfl6BVuD6U5ko5zP\n5+khE+B24E7gpWTl8afIRh8vB26tWfdbwDER8eF8X/PJyunNqvf5kvqMAVKSni7R3NXGCSCldGdE\n/CnwQbJzDB8nu3r78Zr1h4E9ya7yng3cBpyeUno8Io4APkZ2DuR64G/IRh6PAl6Rf9YP8tfbdiKl\nkYh4LXA+cDpZAL0RuLlOnysX1xyf7/OLVcvG+7mb/V4klVyk5O8CSZIkNc/b+EiSJKklBkhJkiS1\nxAApSZKklhggJUmS1BIDpCRJklpigJQkSVJLDJCSJElqiQFSkiRJLfn/gOpapWWDQRYAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x9722668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = range(1,101)\n",
    "y = lst_result\n",
    "Fig, Ax = subplots(figsize=(10, 6), dpi=72)\n",
    "Ax.plot(x,y,'k',lw=2)\n",
    "#Ax.plot(x,[0]*100,'m--')\n",
    "Ax.plot([y.index(min(list(np.abs(y)))), y.index(min(list(np.abs(y))))], [0, min(list(np.abs(y)))],'mo-')\n",
    "Ax.plot(y.index(min(list(np.abs(y)))),min(list(np.abs(y))),'mo')\n",
    "Ax.spines['top'].set_visible(False)\n",
    "Ax.spines['right'].set_visible(False)\n",
    "Ax.yaxis.set_ticks_position('left')\n",
    "Ax.xaxis.set_ticks_position('bottom')\n",
    "Ax.set_title('The function of F(l)')\n",
    "plt.xlabel('The scale of l')\n",
    "plt.ylabel('The values of F(l)')\n",
    "#plt.yticks((-600000,-400000, -200000, 0, 200000, 400000, 600000, 800000, 1000000), ('-600k', '-400k', '-200k', '0', '200k', '400k', '600k', '800k', '1000k'))\n",
    "plt.savefig(r'D:\\result2')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 580,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 580,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_1 = 0.7\n",
    "w_2 = 0.1\n",
    "w_3 = 0.1\n",
    "F_w = np.array([w_1,w_2,w_3,1-w_1-w_2-w_3])\n",
    "F_l_p = (w_1*np.e*mu_L+ w_2*10*average_N_t + w_3*100*np.e*r_k.mean()+ (1-w_1-w_2-w_3)*100*D_L).subs(p_0, 0.01)\n",
    "F_l_res = expand(F_l_p)\n",
    "lst_F_l_res = []\n",
    "for i in range(1,101):\n",
    "    lst_F_l_res.append(F_l_res.subs(l,i))\n",
    "len(lst_F_l_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 581,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import copy\n",
    "F_y_1 = copy.deepcopy(lst_F_l_res)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 578,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(49, 1.53949760319340)"
      ]
     },
     "execution_count": 578,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_y_0.index(min(list(np.abs(F_y_0)))), min(list(np.abs(F_y_0)))  # w_4=0.7\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 582,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(29, 4.16828934755797)"
      ]
     },
     "execution_count": 582,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_y_1.index(min(list(np.abs(F_y_1)))), min(list(np.abs(F_y_1)))    # w_1=0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 570,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(38, 0.989333258474363)"
      ]
     },
     "execution_count": 570,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_y_2.index(min(list(np.abs(F_y_2)))), min(list(np.abs(F_y_2)))   #w_2=0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 584,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn4AAAFwCAYAAAAmMnaNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XmWd///X506abumSFuiqlAotiyyllEWspI4Mg+Bg\nf4qdEQdkcHCcEawyKjMoFvDLoGwj+FVnBFkUcL4uCAp1A8Gx0FB2ZMpW2kILFOidNl2SNsl9/f7I\nnZK2SZvlvpM75PWcx3mcO9d1nXNfJ2dK3l7nXOdESglJkiS9/WX6ugOSJEnqHQY/SZKkAcLgJ0mS\nNEAY/CRJkgYIg58kSdIAYfCTJEkaIIoa/CJidETcHBHrImJDRNyfLz82Ip6MiIaIeCQiZrTZplfr\nJEmSBopij/j9APg48H3gc8DzETEY+BkwHJgPjAN+Gi2G9GKdo52SJGlAKS/WjiNiKvBh4EfAvwG5\nlNIPImIusBfwxZTS9yJiAvBVYA4wqhfrqoF7i3X8kiRJpaZowQ84ML8+EtgENEfEt4A1+fLVO6yn\nAiN7sW6fLhyLJElSv1fMy52D8+thwMeARcCX2DlsRn7d3rvjilknSZI0oBRzxO/F/Pp/Ukq/iIhx\nwPt5K3i9I7+e1Kb9qF6u205EpK997Wvbfq6urqa6urr9o5MkSepbXR7MipTaGxQrjIh4AhgPfAU4\nCzgcmAHcA2wGLs/XNQD7AhXAyt6qSzscfETsWCRJklSquhz8ij2z9W+BZcA1wGjg9JTS08CpwEbg\nP4DXgFNTiy29WVfkY5ckSSopRR3x628c8ZMkSf1IyY34SZIkqUQY/CRJkgaIYs7qlSRJJSrCp5v1\nJ4W6Fc3gJ0nSAOV97f1DIUO6l3olSZIGCIOfJEnSAGHwkyRJGiAMfpIkSQOEwU+SJL1t/elPf6Ks\nrIxMJsPcuXOL8h0rVqzglFNOobKyktGjRzNv3jzWrFnTYfsFCxaQyWTaXV566aWi9LGVs3olSdLb\n0oYNGzj99NMpKyujqampKI+wyeVynHTSSSxdupQTTjiBhoYGfvKTn/Dyyy/zwAMPtLvNMcccw/z5\n87f9vHTpUn7zm98watQo9txzz4L3sS1H/CRJUsm44ooryGQynHPOOQBcffXVZDIZzj//fAAuueQS\nMpkMX/ziF3e7r3PPPZeGhgbOPvvs3batra1l/vz5HS5Llixpd7s777yTpUuXcsghh7Bw4ULuuece\n9t57bxYvXsz999/f7jYnnHACV1111bZl2LBhAJx11lkMHTp0t33tCUf8JElSyaiurgbgwQcfBGDR\nokUA20bPWn+eM2fOLvdz++23c9NNN3HXXXdRU1Oz2+9dv34911xzzU6jgiklIoLDDz+cWbNm7bTd\nY489BsDMmTMByGQyHHbYYaxcuZInnniC4447bpffu2LFCu644w7Ky8s599xzd9vPnjL4SZKk7TTl\nmlhVt6og+5o8cjLlmc7HjcMPP5zKykqefPJJNm/ezKJFi5g2bRoPP/wwDQ0NLF68mLKyMmbPnt3h\nPl577TXOPvtsPvOZz3DiiSd2KvhNmTKFXC7X6X62ar2Xr7KyclvZ8OHDt/Vjd6699lpyuRxz587l\nne98Z5e/v6sMfpIkaTur6laxz7f2Kci+ln9uOVNGT+l0+0wmw+zZs1m4cCG33nora9as4eKLL+bT\nn/40N954I3V1dcycOZMRI0Z0uI/f/OY3rF27lmeeeYaTTz6Z5557DmgZRfzUpz7Fddddt9M2tbW1\nXHTRRR3u87TTTmt3xG/cuHEAbNy4cVtZ6+fx48fv8lg3btzI9ddfT0Rsd89fMRn8JEnSdiaPnMzy\nzy0v2L666rjjjmPhwoVcddVVVFRU8IlPfIILLriAK6+8clt9Z9x3333AW6+mW7NmDffee2+7bdte\n6t3xVXa7utQ7Y8YMAB566CEAmpubefTRRwE49NBDAXj11VdZv349Y8eO3W7yxg033EBdXR1HHHEE\nxx57bKeOqcdSSi75peXXIUnS218p/82rqalJEZEiIh1zzDEppZROOeWUbWW//OUvu7S/BQsWpIhI\nc+fOLXhfc7lcOuCAA1JEpOOPPz69733vSxGRjj766G1tzjjjjBQRaf78+dttt++++6aISLfccssu\nv2MX56rLWcdZvZIkqaTMnDmTyspKImLbSFjrenf39/W2iODuu+/m5JNP5sEHH+Txxx/nox/9KLff\nfvt2bVqXVnfddRfLli1j4sSJzJs3r/f6m3YYzhzIIiL5+5AkDQTtXdJUadrFuerygwm9x0+SJPVL\nF198MdlsdruyiODCCy+kqqqqj3pV2hzxa8MRP0nSQPF2GPHbZ599WLly5XZlEcHy5ct75dEovaWQ\nI34GvzYMfpKkgeLtEPwGikIGPyd3SJIkDRAGP0mSpAHC4CdJkjRAGPwkSZIGCIOfJEnSAGHwkyRJ\nGiAMfipJKSVeyL7Q192QJPVTt912G7NmzaKiooJMJsOZZ55ZtO9asWIFp5xyCpWVlYwePZp58+ax\nZs2aDtsvWLCATCbT7vLSSy8VrZ/QC8EvIoZExLMRkYuIa/Nlx0bEkxHREBGPRMSMNu17tU6l6VfP\n/Yp3f+fdbGna0tddkST1Q0899RQVFRXst99+ANu9J7eQcrkcJ510Er/85S+ZPXs2M2bM4Cc/+Qlz\n587tcJtjjjmG+fPnb1tOOOEEAEaNGsWee+5ZlH626o0RvwuBSfnPKSIGAz8DhgPzgXHAT6PFkF6s\nc7SzhC16eRFbmrfwYu2Lfd0VSVIvuuKKK8hkMpxzzjkAXH311WQyGc4//3wALrnkEjKZDF/84hd3\nuZ9LL72URYsWcfzxx3fqe2tra7cLYzsuS5YsaXe7O++8k6VLl3LIIYewcOFC7rnnHvbee28WL17M\n/fff3+42J5xwAlddddW2ZdiwYQCcddZZDB06tFP97a6ivqs3Ig6hJWx9Fbg8X/xBYC/giyml70XE\nhHz9HGBUL9ZVA/cW8/jVfTWrawB4du2zHLDnAX3cG0lSb6murgbgwQcfBGDRokUAPPDAA9v9PGfO\nnIJ+7/r167nmmmt2GhlMKRERHH744cyaNWun7R577DEAZs6cCUAmk+Gwww5j5cqVPPHEExx33HG7\n/N4VK1Zwxx13UF5ezrnnnlugo+lY0YJffkTtOuDbwMNtqqbk16t3WE8FRvZi3T6dOxL1tuZcMw+/\n0vL/Ms+tfa6PeyNJA0+uKceWVYW51Wbw5MFkyjt/ke3www+nsrKSJ598ks2bN7No0SKmTZvGww8/\nTENDA4sXL6asrIzZs2cXpH+tpkyZQi6X6/J2rffyVVZWbisbPnw4AK+99tput7/22mvJ5XLMnTu3\nV94vXMwRvzOBvYEfAofky0YDFTu0a43W7b2Erph1KlH/+8b/snHrRt77zvca/CSpD2xZtYWafWoK\nsq+jlh/F0Cmdv3yZyWSYPXs2Cxcu5NZbb2XNmjVcfPHFfPrTn+bGG2+krq6OmTNnMmLEiIL0r1Vt\nbS0XXXRRh/WnnXZauyN+48aNA2Djxo3bylo/jx8/fpffuXHjRq6//noigvnz53en211WzOA3GdgT\neKJN2SeA1pu23pFft97/9yItl2V7s24nCxYs2Pa5urp625Czek/N6homjpjI+6e8n3tXeDVeknrb\n4MmDOWr5UQXbV1cdd9xxLFy4kKuuuoqKigo+8YlPcMEFF3DllVduqy+0tpd6U9p+3GhXl3pnzGiZ\nL/rQQw8B0NzczKOPPgrAoYceCsCrr77K+vXrGTt27HaTN2644Qbq6uo44ogjOPbYYwt+TO1KKRVl\nAQ4A/r/8ciGQA+4C3gO8Rkvw+gwtl16X0TISN7g369rpc1Lf+9Qdn0pzfzw33fLkLWmvy/fq6+5I\n0ttSKf/Nq6mpSRGRIiIdc8wxKaWUTjnllG1lv/zlL3e7j9tvvz2dccYZ6YADDkgRkfbdd990xhln\npOuuu66gfc3lctu+4/jjj0/ve9/7UkSko48+elubM844I0VEmj9//nbb7bvvviki0i233LLL79jF\nuepyPivazNaU0tKU0s9TSj8HWqe1LEspPQCcCmwE/iMfyk7NH8CW3qwr1rGrZ2pW13D05KOZNnYa\nr296nXUN6/q6S5KkXjRz5kwqKyuJiG0jYa3rzt7f98QTT3DzzTfz7LPPEhG8+OKL/PCHP9w2OaRQ\nIoK7776bk08+mQcffJDHH3+cj370o9x+++3btWldWt11110sW7aMiRMnMm/evIL2aZf9Nf+8JSLM\ng31s49aNjLpsFPeefi8zJsxg1GWjqPlUDUdOOrKvuyZJbyvtXdJUadrFueryvIWiPs5F6qrW2bwz\nJ86ksqKS8ZXjefbNZw1+kqSdXHzxxWSz2e3KIoILL7yQqqqqPupVaTP4qaTUrKrh3Xu9m8qKlmnx\n08ZOc2avJKldN9xwAytXrtyuLCL4/Oc/b/DrgMFPJaVmdQ1HTXprJtn0sdN5du2zfdgjSVKpWr58\neV93od/xtWUqKTsGP0f8JEkqHIOfSsaqulW8suEVjpq8/Yjf89nnyaWuP01dkiRtz+CnkrF41WIq\nKyo5YI+33s07bew0NjduZnXd6l1sKUmSOsPgp5JRs6qGWRNnUZYp21Y2tWoqZVHm5V5JkgrA4KeS\nseP9fQCDygYxtWqqEzwkSSoAg59KQlOuiUdefWS7+/taOcFDkqTCMPipJPz59T+zuXHzTiN+4CNd\nJEldd8EFF3DooYcyatQoRo0axXHHHVfw17W1WrFiBaeccgqVlZWMHj2aefPmsWbNmg7bL1iwgEwm\n0+7y0ksvFaWPrXyOn0pCzaoa3jHyHUwYMWGnumljp/GLZ3/RB72SJPVXt956K6NGjeJjH/sYNTU1\n/M///A8f/OAHeeaZZ5gwYee/Nd2Vy+U46aSTWLp0KSeccAINDQ385Cc/4eWXX+aBBx5od5tjjjmG\n+fPnb/t56dKl/OY3v2HUqFHsueeeBetbewx+Kgk1q2vavcwLMH2P6axYt4ItTVsYXD64l3smSepN\nV1xxBV/60pf453/+Z6699lquvvpqzjvvPL70pS9x2WWXcckll/C1r32N8847j8svv7zD/dx2220c\nffTRAGzatIlx48axYcMGampq+PCHP7xT+9raWi666KIO93faaacxa9asncrvvPNOli5dyiGHHMLC\nhQvJ5XK8613vYvHixdx///0cd9xxO21zwgkncMIJJ2z7+SMf+QgAZ511FkOHDu34l1MABj+VhJrV\nNZx52Jnt1k0bO41cyrGsdhkH7nlgL/dMktSbqqurAXjwwQcBtl2ebR09a/15zpw5u9xPa+iDllG5\nxsZGIoLJkye32379+vVcc801RMR25SklIoLDDz+83eD32GOPATBz5kwAMpkMhx12GCtXruSJJ55o\nN/i1tWLFCu644w7Ky8s599xzd9m2EAx+6nN1W+pY+sbSdu/vA5hQOYHKikqeW/ucwU+Sesm6dS3L\njkaPblk6076jtrty+OGHU1lZyZNPPsnmzZtZtGgR06ZN4+GHH6ahoYHFixdTVlbG7NmzO7W/pqYm\nPvnJT9LY2Mi8efM44ogj2m03ZcoUcrmuvyyg9V6+ysrKbWXDhw8H4LXXXtvt9tdeey25XI65c+fy\nzne+s8vf31UGP/W5JauXkIkMMyfObLc+Ipg2dhrPvukED0nqLf/xH9Delc+vfQ0WLOhc+47a7kom\nk2H27NksXLiQW2+9lTVr1nDxxRfz6U9/mhtvvJG6ujpmzpzJiBEjdruvzZs3c+qpp7Jw4UJOPvlk\nbr755g7bdvdS77hx4wDYuHHjtrLWz+PHj99l/zZu3Mj1119PRGx3z18xGfzU52pW13DwuIMZNmhY\nh22mj53uI10kqRfNnw+f/OTO5R2N4LXXvqujfa2OO+44Fi5cyFVXXUVFRQWf+MQnuOCCC7jyyiu3\n1e9ONpvl5JNPZvHixZxxxhlcd911lJWVddi+7aXelNJ2dbu61DtjxgwAHnroIQCam5t59NFHATj0\n0EMBePXVV1m/fj1jx47dbvLGDTfcQF1dHUcccQTHHnvsbo+pEAx+6nPtPbh5R9PGTuP3L/6+l3ok\nSerqZdruXNbtSOt9fs888wxHH300Q4cO5T3veQ933nnndvW78qEPfYjFixdTVVXFqFGjOO+88wA4\n8cQTt5tY0aq7l3pPOeUU9t9/f55++mn+8i//ki1btrBq1SqOOuqobQH1X//1X7n55pv53Oc+x9VX\nXw203Dt4zTXXAPTaaB/4HD/1sZQSNas6F/wc8ZOkgWHmzJlUVlYSEdtGwlrXnb2/75VXXiEiWLdu\nHddccw3XXHMN1157LTU1NQXta0Rw9913c/LJJ/Pggw/y+OOP89GPfpTbb799uzatS6u77rqLZcuW\nMXHiRObNm1fQPu2yvzsOZw5kEZH8ffSuletWMuVbU3j6n57e5cSNR155hCO+fwTZL2WpGlrViz2U\npLen9i5pqjTt4lxFe4W74qVe9ama1TWMHDyS/ffYf5ft9hu7HwDPrX2uw+f9SZIGlosvvphsNrtd\nWURw4YUXUlXlIEF7DH7qUzWrajhy0pFkYtd3HYwcPJIJlRMMfpKkbW644QZWrly5XVlE8PnPf97g\n1wGDn/pUzeoaqqdUd6rttLHTfGevJGmb5cuX93UX+h0nd6jPNDY38sirj+x2YkcrH+kiSVLPGPzU\nZ556/Skamho6fenWET9JknrG4Kc+U7Oqhimjp7DX8L061X76HtN5fu3z5FLXn7MkSZIMfupDnXlw\nc1vTxk6jvqmeVXWritgrSZLevpzcoT5Ts7qGsw8/u9Pt9xm9D+WZcp5b+xzvHFX8F1lL0ttd2wcK\na2Aw+KlPrGtYxzNvPtOlR7MMKhvE1KqpPPvms3xg6geK2DtJevvz4c0Dk5d61SeWrF5CeaacGeNn\ndGk7X90mSVL3FS34RcR+EfGHiHgzIuoi4rcRMTVfd2xEPBkRDRHxSETMaLNdr9apb9SsruHQcYcy\ndNDQLm03fex0nssa/CRJ6o5ijvhNzK8vBG4APgBcFxGDgZ8Bw4H5wDjgp9FiSC/WOdrZh7o6saPV\ntLHTePZNH+kiSVJ3FPMevwdSSnNaf4iITwAHAicCewFfTCl9LyImAF8F5gCjerGuGri3iMevDqSU\nqFlVw6l/eWqXt50+djor1q1gS9MWBpcPLkLvJEl6+yraqFdKqbH1c0QcAVQBfwT2yRev3mE9tZfr\nWsvVy1asW8Ebm9/o9ohfIvFC9oUi9EySpLe3os/qjYj9gTuB5cA5wGk7Nsmv25teVMy6di1YsGDb\n5+rqaqqrq3fVXN1Qs7qG0UNGs9/Y/bq87fjK8YyoGMFza5/joL0OKkLvJEl6+ypq8IuIA2m5nLoZ\neH9KaU1EvJivfkd+PSm/fpGWy7K9WbeTtsFPxbF41WKOnHQkmW7cZhkRvrpNkqRuKlrwi4h3AH8A\nxgBXA8dExNHAL4DXgc9ExEbgLFpGA+8DKnq5Tn2gZnUNx089vtvbT99juo90kSSpG4o5s/VdwJ75\n7/h34Fbg1pTSFuBUYCPwH8BrwKmpRa/WFfHY1YGtzVt57NXHunV/X6tpYxzxkySpO8L885aIMA8W\n2ZLVSzjyuiN5/V9eZ8/he3ZrH7c9dRvn/vpc3vjiGwXunSRJ/UqX37nns+zUq2pW1zC1amq3Qx+0\nXOp9c/ObZOuzBeyZJElvfwY/9aruPri5rf3GtMwG9j4/SZK6xuCnXlWzqufBb8TgEUwcMdHgJ0lS\nFxn81Guy9Vmezz7P0ZOP7vG+fHWbJEldZ/BTr3lo9UNUlFVw2PjDeryv6WOn81zWET9JkrrC4Kde\n89P//SnHTD6mIO/YdcRPkqSuK/or2ySAVXWruPmJm/nZx35WkP1NHzud57PPk0u5br0BRJKkgci/\nmOoVVz5wJdP3mM5J004qyP6mjZ1GQ1MDq+pWFWR/kiQNBAY/Fd0bm97gvx79L/71vf9asNG5KaOn\nUJ4p93KvJEldYPBT0V1Tcw3jK8fzsYM+VrB9DiobxLuq3uUjXSRJ6gKDn4qqbksd1z50LV8+9suU\nZwp7S+m0sb6zV5KkrjD4qai+9/D3GDZoGGccekbB9z197HRH/CRJ6gKDn4qmvrGeqx68in95z78U\n5BEuO3LET5KkrjH4qWhuePwGGnONnD3z7KLsf/oe01m5biUNTQ1F2b8kSW83Bj8VRWNzI99c9E0+\nd9TnqKyoLMp3TBs7jUTihewLRdm/JElvNwY/FcVtf76NtfVr+eyRny3ad4wbPo6Rg0d6n58kSZ1k\n8FPB5VKOy/50Gf848x8ZM3RM0b4nIpg2dprBT5KkTjL4qeDueOYOltUu4wvHfKHo3+UED0mSOs/g\np4JKKXHpny7l7w/7eyaMmFD07/ORLpIkdZ7BTwX1+xd/z2OvPsYXj/1ir3zftLHTfG2bJEmdZPBT\nQV36p0v524P/lqlVU3u0n+ZNzaz707rdtps+djpr69eydvPaHn2fJEkDgcFPBfPAyw9w34r7OP/Y\n83u0n9r7allyyBIen/04j335JXK51GHb/cbuB+DlXkmSOsHgp4L59z/9Ox/e/8MctNdB3dq+aWMT\nz332OZ54/xNUfaCK/X50AH99+R6cfeS6DsNfZUUlk0ZMMvhJktQJ5X3dAb09PLnmSX713K+o+VRN\nt7a/86oNrL1yBfuXb+SQ3x7CmA+0PAbm2lc28LdfGsnGg9dzyxOjKCuPnbZ1Zq8kSZ3jiJ8K4rI/\nXcZf7PMXHDnpyC5tt+61Jk47eB0fPq+SpyZPYNafZ20LfQAf/uIIfvHden71v5XMnVbH1obcTvtw\nZq8kSZ1j8FOPvZB9gf9++r/5t9n/1qXt7rhyAwe+o4l7lw7hZ9/YyFU1e1A+YudB6BP+sZKFP9rC\nH1cM44NTN9BUv3348yHOkiR1jsFPPfbNRd9k1sRZzJkyp1PtmzY08ZmZa5n7L5VUT2/gf18qZ+6X\nRuxym9mnDed3v2jk4E1reXruUzRvbt5WN23sNJ7PPk8u7TwaKEmS3mLwU4+srlvNjY/fyL/N/jci\ndr7/bke199Sy5OAlTF31JrdfsZFb/zyaqomdu9V01l8P49JHJ7D5mc08eeKTNNU1AXDgngfS0NTA\nR/7fR7jlyVtY17D7x8BIkjQQRUodPyrj7SQijgW+C0wDngY+lVJ6bIc2aaD8PnqiKddEzaoafv3C\nr/n5Mz8nExme+McnyET7/zti87ObWfurtbz5yzdZ/8f1TPyniUy9bCrlld2bW9TwcgNPfOAJykeV\nc8ivD2HQmEEsfH4htzx1C7967ldsatzEnClzmLv/XE7Z/xQmjpjYk8OVJKlU7X7EZccNBkLQiYgh\nwApgE3A58BVgC7BfSm9dHzT4deyl9S/xmxd+w6+X/Zp7XryHui11HDHxCP5q37/izMPOZJ+qfba1\n3bI5x2//axO/+FETf3yqgmu3PkLVvoMZ+6Gx7PU3ezHyyJE97s/WNVt54i+fgASH/u5QKsZVtJQ3\nb+X+Ffdz+zO384tnfsGrG1/lqElH8eH9P8zc/ecyfY/pPf5uSZJKhMGvPRExF/gZ8MWU0pURcRHw\nVeADKaV727Qz+OXVN9Zz/8r7t4W9Z958hvGV4znhXSfwV/v+FR+Y+gH2GLbHtvaNtY386OJN/PTn\nwR9fHs7mVMahIzfxV+9t5JyLKhg/c1inLgV3RWO2kbvnPMNFz+/Nf/+hgv2OGrJdfS7lWLJ6Cbc/\nczu3P3M7z619jgP2OIBTpp/CtLHTGDtsLGOHjmXM0DGMHdayLs/4hCNJUr9h8GtPRHwBuAL4eErp\nxxFxNvA94B9SSte3aTeggl8ul9i8qZGVbyxn+drnWfb6Ml58/UVWvL6CpauXkt6cwr6ZOew35FAm\nV0xjaPMebN4IR0yq510jt5DbnKN5czMbH9/I+j+t58rMdJrGD+OkkxIf/cJwxu83qOjHkF3dxPEH\n1bNyw2AO2aue4UMSw4bAXx+8kffs30jZ8DLKKsvIDMvwWnqN25cuZfHqZWSbX2dD0zoa0iZSpolM\n5SvE8LUMHzycEUNHMGrYKEYNG8Xg3HjK0ygGlZdTUV5BRfkgBg0qZ+SIROXwDIPKB1FRXsHg8sFU\nlFeQaxxCc2M5ZZkyMpkM5WVllEUZQ4fCkMFBWabsrSXK2NqYITWXEQSZTMul8ohgcAUMHpwhIshE\nhiCITNC4Ncg159sRRARBUDEYysvblOVDdmNj0NzcUt667yAYNAgGDXqrrFVjI+Ry2/93JAjKy/P7\n3yG8t7R/a/+tysuhfNDO/z1qbISU27m8Zf/tt8+1M2dnUHlQVrZzeUfty8vpVPvW4+ts+87uv73f\nTzH6v7v2TU3tty8rs/1Abr/jv+u+7o/tO99+5OghXQ5+A3V4o7BDT0WyZXOON1c2Ub65iSGNjTRl\nm2iqbaKxtuXz7YuH8NDKIWzZClsao2XdFJw64U3eW7mO3JYcuYaWpbmhmWuy7+Tu+vE0kaGJoIkM\nUMHnCf6aSg7N/1+rq5jGL5nIPeQYEs0MyeQYmsnxjxPXUjVpPWXDysgMzzBixgimXDiFO2aPIjOo\nd+cLjZlUzn0vDOOS0+p4/XXYuDnYsDmoW9nI+jfW07ypmeaNzdvWT607jt/lPr7TfubzHKfwyk7l\nV7MfdzKp3fbv6WL7U3iFHDkaaex0+67u3/a2t73tbT9w2ndnrGqgjPh9GPg58OWU0uURcTEt9/n9\nRUrpD23apa997Wvbtquurqa6urro/WtY2cC6P67jv/4LbntkBLVby9nQXMbmfC7/LM/zEVYTg4NB\nYwZRPqacQVWDuKNuTxbXjmDIoERFBQypgMGD4S/2X8/EcSt5qeElVtSvYNmmZazespps3T5kcvsw\nduQI9hw1inFVY5g4Zg+OPHAMB79rL8oGlxEVQWZQhszQDI1lGSoqyxg0uF/k5E7b2pCjsSHRvBWa\nGhPNTYkh5TB4UCLlEjRDam75/PrrsH49kCDlIJcSKQd7jUlUjYSUEiS2LavXBNl1QUqJbf+yEkzc\nM0fV6GYArr2wAAAZqklEQVSam5vJkSOXy9Gca+blVzO8mW9Park8nUhM2LORsaNb2iZa+pJSYtVr\n5axdV0bbf7eJxIS9tjK2qolEgtxb5atfq2Btbdm2n1t7NXFcvv0Or8J7+bVB1K57a6S29XsmjtvC\nmKpGdrT6tcFka3ce2Z2wVwNj22m/6rUKsuvaaT+u/farXxu8c/uUmDBuS4fta9vdf9+1b/f3s4vj\nrV1XUbz2rw4hu37n/ky0/YBtv3b9zuM//an/A739Z84/1Eu97YmIwcBKYDNvTe5oAPZte223Ny71\n5nKJNY9sJj26nnX/s471f1zPlpe3UD62nCf2nczKUaOYMCnYY3yGMROCPSeXMWX/MsZOKads6M5j\nv43NjTz1+lMsWb2Eh1Y/xJJXlvD0G0+TUmL/PfZn1qRZzJrYshy010FUVlQW9fgkSVKvMfh1JCJm\nA/8XmA78mZb7+x7doU3Bg18ul1h022Z+/9OtLKrJ8Mhrw5iZslwy+UVGvW8Uo2aPYvTs0Qw7YBiR\n6dz5a84189P//SnfXvJtlqxewpbmLew9au9tIe/ISUdy+ITDGTm457NnJUlSyTL49UQxgt9F/986\nFtw+mncOqufIfbYw+zg44fQhTDt2cJdnuW5t3sqPnvwRl/3pMl6ue5mzZpzFifueyKxJs9hr+F4F\n7bckSSp5Br+eKHTwa97czPuq1nPcB4JL76rq9n7qG+u57tHruPyBy1nXsI5/nvXPfP6Yzxv2JEka\n2JzVW0peve5VLqtcwaxbju7W9nVb6vjuku9y1eKraMo1Mf+o+Xz2yM9SNbT7IVKSJA1cBr8iyW3J\n8dI3X+Kd572DIaO79mteu3kt36r5Ftc+dC1Dy4fy5WO/zNkzz3ZihiRJ6hGDX5G8esOr5DblmPTZ\nnZ/R05GGpga+eu9X+e7D32Wv4Xtx2V9cxhmHncGQ8iG731iSJGk3DH5FkGvM8dJlLzHpc5MoH9n5\nX/Gl/3MpNz1xE9896bv87cF/6+vDJElSQZksiuC1H66haW0Tk8+d3OltlmWX8c1F3+QHp/yAjx+8\n85slJEmSespZvW0UYlZv45Ych43azPlzN/B3t03o9HYfuu1D1G2p474z7uvyY14kSdKA5Kzevnbd\nFzbw4pZK3n/hzq9R6sivnvsVC59fyKOfftTQJ0mSisYRvzZ6OuLX3JTYb3g9c969hesf6dwjVxqa\nGnj3d97NSfudxLdO/Fa3v1uSJA04jvj1pZvOr2PV1hF89fqd36nbkSseuIINWzdw0ZyLitgzSZIk\nR/y205MRv1wuccDwzczar5EfPTm6U9usXLeSA/7vAXznpO/wycM+2a3vlSRJA5Yjfn1l+Y+z7LEl\nceF/dv4hy+f99jwOHX8opx96ehF7JkmS1MLgVwApJdb9xwp+cFYl047Zo1Pb/G7Z7/j50p+z5B+W\nkIlMkXsoSZJk8CuI2t/WsuHRDRz44wM71X5r81bOWXgOn575aWZOnFnk3kmSJLVwqKmHUkqsuGQF\n404bx9CpQzu1zbcWf4s3Nr/B19//9SL3TpIk6S2O+PXQuvvXUfdAHftfv3+n2q+uW83Ff7yYK46/\ngrHDxha5d5IkSW9xVm8b3ZnVe+/spxg3KcNBPz6oU+0//rOP8+zaZ3noUw9Rlun8Y18kSZJ24Kze\n3vTb/9zIyX86iKX31neq/f0r7ufHf/4xD5z1gKFPkiT1Okf82ujqiN9796xj8CC455WRu23blGti\nxn/OYNbEWfzglB/0pJuSJEngiF/vue/mjSx6cyT3/3BTp9p/Z8l3eHn9y/z+735f5J5JkiS1zxG/\nNroy4vf+Cetpag7++PruR/vWbFzDtG9P45I5l3DuUef2tJuSJEngiF/vWPzTzfzhtVH87vsbO9X+\n/HvOZ+9Re/NPs/6pyD2TJEnqmMGvG7b85nXm71vBBz41cbdtH37lYW56/Cbu++R9lGf8dUuSpL5j\nEumG9/3n3ryntqlTbR98+UEO2PMA3rf3+4rcK0mSpF3zzR3dEJlg0NhBnWqbrc8ydqgPapYkSX3P\n4Fdk2fosY4aO6etuSJIkGfyKrbah1uAnSZJKgsGvyBzxkyRJpaIowS8ivh0RKyKiPiKejYiPt6kb\nERG3RcSmiHg1Is7rq7reYPCTJEmlolgjfkcANwBfAEYDN0XEPvm6rwPzgG8ADwKXR8ScPqorumx9\nlqohVb31dZIkSR0qyps7ImJQSqkx//lK4PPAB1NKv46IdcDLKaWD82FwGfDDlNIZvVj3o5TS6e30\nu0vv6u2MvS7fi2tPvJZ5755X0P1KkqQBr8tv7ijKiF+b0DcImANsAh6JiDHASGB1vmnrempEVPVi\nXevoY1GllLzUK0mSSka3g19ErIqIXDvL6fn6cuBHwCHAP6SU3mhvN7v6il6uK7iNWzfSnJoNfpIk\nqST05M0ds4H2nmL8Wn6k78fAh2kJfT8GSCllI6IOmJxvOym/fjGlVBsRG3qrrqODWrBgwbbP1dXV\nVFdXd9R0t7L1WQCDnyRJKgndDn4ppeUd1UXEbcBc4C5gU0T8DbA4pbQCuAk4JyIuBGYACbgxv+mN\nvVy3k7bBr6cMfpIkqZQU6129R9MSsE7KLwk4E1gBfAUYB3wZqAPOTyn9Ib9db9cVVbY+SyYyjBg8\noje+TpIkaZeKEvxSSh1OnkgpbQD+phTqiq31US6Z8DnZkiSp75lIisjXtUmSpFJi8CsiH+UiSZJK\nicGviAx+kiSplBj8iihbn6VqqK9rkyRJpcHgV0TZ+ixjhjjiJ0mSSoPBr4i81CtJkkqJwa+InNUr\nSZJKicGviBzxkyRJpcTgV0QGP0mSVEoMfkXS0NTA5sbNzuqVJEklw+BXJLX1tQCO+EmSpJJh8CuS\nbH0WMPhJkqTSYfArktqGlhG/qiFe6pUkSaXB4Fck2fosIypGMKhsUF93RZIkCTD4FY2va5MkSaXG\n4FckPspFkiSVGoNfkRj8JElSqTH4FUltva9rkyRJpcXgVyTZhixjhhj8JElS6TD4FYmTOyRJUqkx\n+BWJ9/hJkqRSY/ArEoOfJEkqNQa/IjH4SZKkUmPwK4LmXDPrG9Yb/CRJUkkx+BXB+i3rSSSDnyRJ\nKikGvyLI1mcBqBrirF5JklQ6DH5F0Br8HPGTJEmlxOBXBNn6LBVlFQwbNKyvuyJJkrRNUYNfRFwc\nEbmI2NCmbERE3BYRmyLi1Yg4r6/qiqV1Rm9EFPurJEmSOq1owS8iDgL+BWgAUpuqrwPzgG8ADwKX\nR8ScPqorCt/TK0mSSlFRgl9EZIDrgP8E1uxQfQbwdErpYqB19O2TvVx3ZjcPrVOy9VkndkiSpJJT\nrBG/fwbGAV8Btl3vjIgxwEhgdb6odT01Iqp6sW6fHh3dbvjwZkmSVIq6HfwiYlX+/r0dl88BlwKX\nAxOA8pbmMbW93ezqK3q5rmCyDQY/SZJUesp7sO1sYFA75UOA4cD/3aH8mZRSRX6ix+R82aT8+sWU\nUm1v1nV0UAsWLNj2ubq6murq6o6adihbn2W/Mft1eTtJkqRi6nbwSyktb688IoYCp9IyoSOA7wKV\nwMfzTW4EzomIC4EZ+XY39lHdTtoGv+5ycockSSpFPRnxa1dKqR74WevPEXEFMCSl9It80Vdouf/v\ny0AdcH5K6Q99VFcU3uMnSZJKUcGD345SSvvs8PMG4G86aNurdcXirF5JklSKfHNHgaWUHPGTJEkl\nyeBXYJsaN9GYazT4SZKkkmPwK7BsfRbA4CdJkkqOwa/AautrAYOfJEkqPQa/AsvWZwmCUUNG9XVX\nJEmStmPwK7BsfZbRQ0aTCX+1kiSptJhOCswZvZIkqVQZ/ArM4CdJkkqVwa/ADH6SJKlUGfwKrLah\nlqqhvrVDkiSVHoNfgWXrs4wZ4oifJEkqPQa/AvNSryRJKlUGvwIz+EmSpFJl8Cswg58kSSpVBr8C\nq22oNfhJkqSSZPAroK3NW9m4daOzeiVJUkky+BVQbX0tgCN+kiSpJBn8CihbnwUMfpIkqTQZ/Aqo\nNfhVDfFSryRJKj0GvwLK1mcZPmg4g8sH93VXJEmSdmLwKyBf1yZJkkqZwa+AfIafJEkqZQa/AjL4\nSZKkUmbwKyCDnyRJKmUGvwLK1mcZM8TgJ0mSSpPBr4Ac8ZMkSaXM4FdAzuqVJEmlzOBXQI74SZKk\nUla04BcR/xgRyyKiISJeiIj35stHRMRtEbEpIl6NiPPabNOrdYVm8JMkSaWsKMEvIv4a+A7wEvBP\nwK3AoHz114F5wDeAB4HLI2JOH9UVTC7lqK2vNfhJkqSSFSmlwu804o/ADGASsDWl1NCmbh3wckrp\n4IjYB1gG/DCldEYv1v0opXR6O/1O3f191NbXMuabY3js049x2PjDurUPSZKkLoiublCsS70HAo3A\nUmBTRCyKiEkRMQYYCazOt2tdT42Iql6s26cAx7id2oZaAKqGOLlDkiSVpm4Hv4hYFRG5dpYzgMHA\naOAa4ALgGOAyYMfhtF0l1d6u65FsfRbAS72SJKlklfdg29m8dd9eW68B5wEHAVfREi4vBaamlGoj\nYgMwOd92Un79Ym/XdXRQCxYs2Pa5urqa6urqjppuJ1ufpTxTTmVFZafaS5Ik9bZuB7+U0vKO6iLi\nRuAKWgJf6yjbH/PrG4FzIuJCWu4DTPmyvqjbSdvg1xWtM3ojijaoKEmS1CM9GfHblWuAdwH/AGwF\nvg9clK/7CjAO+DJQB5yfUvpDH9UVjI9ykSRJpa4os3r7q57M6v36H7/O3c/fzQNnPVDgXkmSJLWr\nZGb1Djg+w0+SJJU6g1+BZBu81CtJkkqbwa9AvMdPkiSVOoNfgRj8JElSqTP4FYjBT5IklTqDX4Fk\n67O+rk2SJJU0g18BpJSc1StJkkqewa8A6pvq2dK8xeAnSZJKmsGvALL1WQCDnyRJKmkGvwIw+EmS\npP7A4FcArcFv9JDRfdwTSZKkjhn8CqC2vpZRg0dRlinr665IkiR1yOBXAD7DT5Ik9QcGvwIw+EmS\npP7A4FcABj9JktQfGPwKIFufpWqob+2QJEmlzeBXANmGLGOGOOInSZJKm8GvAHxdmyRJ6g8MfgXg\nPX6SJKk/MPgVgMFPkiT1Bwa/AjD4SZKk/sDg10ONzY1s2LrBWb2SJKnkGfx6qLahFsARP0mSVPIM\nfj1UW2/wkyRJ/YPBr4ey9VkAqoZ4qVeSJJU2g18PZeuzDC0fytBBQ/u6K5IkSbtk8OshX9cmSZL6\nC4NfD/koF0mS1F8Y/HqotsHXtUmSpP6hKMEvIo6OiJqI2BgRayPi/0XEyHzdiIi4LSI2RcSrEXFe\nm+16ta4QHPGTJEn9RbFG/L4FzAIuA+4HPgqcm6/7OjAP+AbwIHB5RMzpo7oey9ZnGTPE4CdJkkpf\nsYJfLZAD7gUeb1MGcAbwdErpYqB19O2TvVx3ZvcPbXuO+EmSpP6ivEj7/QfgT/kF4NfAdyJiDDAS\nWJ0vb11PjYiqXqzbp2eH9xZn9UqSpP6i28EvIlYBE9upOhP4SL7uLOBAWkbazgF+uONudvUVvVwH\nwIIFC7Z9rq6uprq6epftHfGTJEn9RU9G/GYDg9opfw34LrA8pXRDREymJfgdn1K6JiI2AJPzbSfl\n1y+mlGp7s66jg2ob/DrDWb2SJKm/6HbwSykt76guIpYCh0XEl4D988XP5tc3AudExIXADCDly/qi\nrkdyKeeInyRJ6jeKdY/fmcC1wFeALcCPgUvydV8BxgFfBuqA81NKf+ijuh7ZsGUDuZQz+EmSpH4h\nUkp93YeSERGpK7+P5bXLmXrNVF4890X2qSrYfBFJkqTO2O3chR355o4eyNZnARzxkyRJ/YLBrwey\n9VnKooyRg0f2dVckSZJ2y+DXA7UNtVQNrSKiyyOtkiRJvc7g1wPO6JUkSf2Jwa8HsvVZqob41g5J\nktQ/GPx6wBE/SZLUnxj8esDgJ0mS+hODXw/4ujZJktSfGPx6wBE/SZLUnxj8esDgJ0mS+hODXw84\nq1eSJPUnBr8ecMRPkiT1Jwa/bqpvrKehqcHgJ0mS+g2DXzfVNtQCGPwkSVK/YfDrpmx9FjD4SZKk\n/sPg102twW/0kNF93BNJkqTOMfh1U7Y+y4iKEQwqG9TXXZEkSeoUg183OaNXkiT1Nwa/bjL4SZKk\n/sbg10219b6nV5Ik9S8Gv25yxE+SJPU3Br9uyjb4ujZJktS/GPy6yRE/SZLU3xj8usngJ0mS+huD\nXzc5uUOSJPU3Br9ucsRPkiT1Nwa/bmjKNbF+y3qqhjq5Q5Ik9R8Gv25Y17AOwBE/SZLUr3Q7+EXE\niRHxVETk8suYNnUjIuK2iNgUEa9GxHmlWNdd2fosYPCTJEn9S09G/IYC9wMvAGmHuq8D84BvAA8C\nl0fEnBKs6xaDnyRJ6o8ipR0zWxd3EHEfMBvYM6WUzZetA15OKR0cEfsAy4AfppTOKJG6H6WUTm/n\nWFJnfh8Ln1/I3P+eS/0F9UREj35/kiRJ3dTlEFLwe/zyl3xHAqvzRa3rqRFRVSJ1+3T/CFtG/KqG\nVhn6JElSv7LL4BcRq9rcw9d22Wm0bFe76Sd1neajXCRJUn9Uvpv62cCgdspf62iDlFI2IuqAyfmi\nSfn1iyml2ojYUAp1HfV/wYIF2z5XV1dTXV29UxuDnyRJ6o92GfxSSss7qouIfYFqYAItI2l/FxHP\np5TuBm4CzomIC4EZtEz+uDG/6Y0lVLeTtsGvI025JiZUTthtO0mSpFLS7ckdEfFJ4Ae8NaM3gPtS\nSu+PiBHA94EPAXXAVSmly/PblUxdO8fUqckdkiRJJaDLt7D1eFbv24nBT5Ik9SN9P6tXkiRJpcng\nJ0mSNEAY/CRJkgYIg58kSdIAYfCTJEkaIAx+kiRJA4TBT5IkaYAw+EmSJA0QBj9JkqQBwuAnSZI0\nQBj8JEmSBgiDnyRJ0gBh8JMkSRogDH6SJEkDhMFPkiRpgDD4SZIkDRAGP0mSpAHC4CdJkjRAGPwk\nSZIGCIOfJEnSAGHwkyRJGiAMfpIkSQOEwU+SJGmAMPhJkiQNEAY/SZKkAcLgJ0mSNEAY/CRJkgYI\ng58kSdIA0e3gFxEnRsRTEZHLL2Pa1H07IlZERH1EPBsRH29TNyIibouITRHxakSc11d1kiRJA0lP\nRvyGAvcDLwBph7ojgBuALwCjgZsiYp983deBecA3gAeByyNiTh/V6W3kvvvu6+suqJs8d/2b56//\n8tz1bxFR3dVtuh38Uko/Tyl9FnilnerZKaWLUkrfBX4ElAHT83VnAE+nlC4GWkffPtnLdWd26WDV\nL/gfsP7Lc9e/ef76L89dv1fd1Q2Kco9fSqkRICIGAXOATcAj+cvBI4HV+aat66kRUdWLda2jj5Ik\nSQPGLoNfRKxqcw9f2+X03e04IsppGe07BPiHlNIb7TXb1S56uU6SJOntLaXU4ULLyNi0dpaRbdrc\nBzQDY9qUDQJ+li//+x32uQ74c/7zVCAH3JT/eX1v1rVzvMnFxcXFxcXFpb8su8px7S3l7EJKaXlH\ndRGxLy3XlifQMpL2dxHxfErpbuBmYC5wF7ApIv4GWJxSWgHcBJwTERcCM/IdvzG/2xt7uW7H43VE\nUJIkvW1FfqSr6xtGfBL4AS1BClrC330ppfdHxHLgnbx1aTUBZ6aUbo6IEcD3gQ8BdcBVKaXL8/vs\n1TpJkqSBpNvBT5IkSf2Lb+4AIuLYiHgyIhoi4pGImNHXfVL7IuKaiFiTn2T0yzblnsN+ICL2i4g/\nRMSbEVEXEb+NiKn5Os9hiYuImvx52xQRj0XEX+XLPXf9REQMyb9YIRcR1+bLPH/9QP7FGG0n2j6W\nL+/S+RvwwS8ihtAyEWU4MB8YB/w0Igb876ZEJeC2Np89h/3LxPz6Qloe8v4B4LqIGIznsD9YBJwD\nXAK8G/iu567fuRCYlP+cPH/9SqLlxRl/k1++1J2/f55YOBHYC/hOSul7wPW0zGau7stOqX0ppc8B\nV+9Q7DnsPx5IKc1JKX0nfy5rgQPxHPYLKaUv0DJp715gCy0P8P8gnrt+ISIOoSUcfK1Nseev/whg\nBXB3Sun/pZR+Rzf+22nwe+thzj7kuf/Ycfa157CfaH24O0BEHAFUAX/Ec9gvRMRo4HVgMS2jD+fh\nuesX8iNA1wHfBh5uUzUlv/b8lb4EnA7U5W95+nu68e/P4LczH+nS/3kOS1xE7A/cCSyn5dLhjufM\nc1iaNgDHA+fS8irOH7fTxnNXms4E9gZ+CEzOl40GKnZo5/krXd8HTgVOAxqA/2ynzW7P3y6f4zdA\nvJhfvyO/nrRDuUpf6/MmPYf9QEQcSMulws3A+1NKayLCf4f9QEqpGbgHuCciTgVmAy/lqz13pW0y\nsCfwRJuyT+DfwH4jpXRp6+f8FZPPA6vyRZ0+fwP+cS75G1tX0vJH6HLgK7Qk6X3TQP/llKCIOImW\nm8r/HXgSuBZ4CPgdnsOSFxHvoOUy0xhaztNKWi5f/AL/HZa0iDgB+BjwAC1/ZC4AXgYOwHNX8iLi\nAFrOFbT8N3QBsBD4P8DP8fyVtPz9mf+HlnNWDnwVGArsBzxGV85fV1/18XZcaPlfrU/ScrPyI8Dh\nfd0nlw7P1R9oee1ec5v16Z7D/rHQcsNx2/OXA5rzdZ7DEl6AI4Cn8n9gsvk/QAd57vrfAhyX/7d3\njeevfyzAeFomVr0BbKJlwOP47py/AT/iJ0mSNFA4uUOSJGmAMPhJkiQNEAY/SZKkAcLgJ0mSNEAY\n/CRJkgYIg58kSdIAYfCTJEkaIAx+kiRJA8T/D5xVemOCPLeMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xb825780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sensitivity analysis\n",
    "fig0, Ax0 = subplots(figsize=(10,6))\n",
    "F_x = range(1,101)\n",
    "p1 = plt.plot(F_x,F_y_0,'g-', lw = 1.25, label='w_4 = 0.7')\n",
    "p2 = plt.plot(F_x,F_y_1,'m-', lw = 1.25,  label='w_1 = 0.7')\n",
    "p3 = plt.plot(F_x,F_y_2,'b--', lw = 1.25,  label='w_2 = 0.7')\n",
    "Ax0.spines['top'].set_visible(False)\n",
    "Ax0.spines['right'].set_visible(False)\n",
    "Ax0.yaxis.set_ticks_position('left')\n",
    "Ax0.xaxis.set_ticks_position('bottom')\n",
    "plt.xlim(0,50)\n",
    "#plt.ylim(0,40)\n",
    "Ax0.legend(loc='upper right')\n",
    "#plt.yticks((-200000,-100000, 0, 100000, 200000, 300000, 400000), ('-200k', '-100k', '0', '100k', '200k', '300k', '400k'))\n",
    "plt.savefig(r'D:\\sensitivity4')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
